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ABSTRACT 

We use the Herschel- ATLAS survey to conduct the first large-scale statistical study 
of the submillimetrc properties of optically selected galaxies. Using ~ 80, 000 r-band 
selected galaxies from 126 deg 2 of the GAMA survey, we stack into submillimetre 
imaging at 250, 350 and 500 /im to gain unprecedented statistics on the dust emission 
from galaxies at z < 0.35. We find that low redshift galaxies account for 5% of the 
cosmic 250 /xm background (4% at 350 /xm; 3% at 500 /Ltm) , of which approximately 
60% comes from 'blue' and 20% from 'red' galaxies (rest-frame g — r). We compare 
the dust properties of different galaxy populations by dividing the sample into bins of 
optical luminosity, stellar mass, colour and redshift. In blue galaxies we find that dust 
temperature and luminosity correlate strongly with stellar mass at a fixed redshift, but 
red galaxies do not follow these correlations and overall have lower luminosities and 
temperatures. We make reasonable assumptions to account for the contaminating flux 
from lensing by red sequence galaxies and conclude that galaxies with different optical 
colours have fundamentally different dust emission properties. Results indicate that 
while blue galaxies are more luminous than red galaxies due to higher temperatures, 
the dust masses of the two samples are relatively similar. Dust mass is shown to 
correlate with stellar mass, although the dust-to-stellar mass ratio is much higher for 
low stellar mass galaxies, consistent with the lowest mass galaxies having the highest 
specific star formation rates. We stack the 250 /j,m-to-NUV luminosity ratio, finding 
results consistent with greater obscuration of star formation at lower stellar mass 
and higher redshift. Submillimetrc luminosities and dust masses of all galaxies are 
shown to evolve strongly with redshift, indicating a fall in the amount of obscured 
star formation in ordinary galaxies over the last four billion years. 

Key words: galaxies: statistics - galaxies: ISM - galaxies: evolution - submillimetre: 
galaxies - submillimetre: diffuse background. 



1 INTRODUCTION 

Dust in galaxies represents only a tiny fraction of the mass 
density of the Univers^E yet from an observational point of 
view it can provide some of the most important indications 
of the history of star formation. This is possible because 
most stars form in dense clouds of gas and dust. The dust 
in these regions is heated as it absorbs ultraviolet (UV) ra- 
diation from the hot massive stars that form within, and re- 
radiates the energy as a modified blackbody (or greybody) 
with a characteristic temperature generally around 20—30 K. 
Measurement of the UV radiation from these massive stars is 
the most direct method of inferring the star-formation rate 
(SFR), although as with all SFR indicators this relies on 
estimating the rate of massive star formation from the UV 
and assuming an initial mass function (IMF) to infer the to- 
tal SFR. Dust poses a problem to this method since the UV 
attenuation must be accounted for, and this will vary from 
one star forming region to another as well as being depen- 
dent on the inclination angle that the galaxy disk is observed 
at. Observing the thermal emission of dust with far-infrared 
(FIR) and submillimetre (submm) telescopes provides a way 
to trace the UV radiation field in all parts of the interstellar 
medium (ISM) of a galaxy, since in general the ISM is op- 
tically thin to FIR wavelengths. Hence one can use the UV 



* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and 
with important participation from NASA, 
f E-mail: ppxnbl@nottingham.ac.uk 

The cosmic dust mass density was estima ted to be 0.0083 pe r 
cent of the baryon density at redshift 0.0 bv lDriver et al. I j2007h 



and FIR in tandem to measure the (massive) SFR in the 
unobscured and obscured phases respectively. 

A further complication arises from the fact that dust 
also exists in the large-scale ISM, in regions not associ- 
ated with star formation. Indeed this component forms the 
bulk of the dust mass in a galaxy, and in spiral galax- 
ies can dominate the bolometric output in the FIR ( Heloul 



unnc fe Ealedl200ll ; ISauvage, Tuffs fc PopescJ : 
Draine et al.l 120071 ). Because this ISM component is heated 
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by the general interstellar radiation field, its FIR emission is 
not necessarily correlated with the SFR. The extent to which 
both the UV and FIR luminosities can be used to trace SFRs 
can only be understood by studying the full spectral energy 
distribution (SED) of galaxies from the UV to the FIR and 
submm. 

We know that the mid/far-IR luminosity density of 
the Universe increases towards higher redshifts, as a re- 
sult of increased star-formation activity and increased 



dust content dBlain et all 1 19991: iFranceschini et all 1200 ll : 
iDunne. Eales fc EdmundsN2003l : iLe Floc'h et al.ll2005l ). Re- 

cent analyses of the submm luminosity function (LF; 



lEales et all 120091 . l2010al : iDve et aLlkoid ; IDunne et al]|201ll 
[hereafter |Plll |) with the Balloo n-borne Large Ape rture 
Submillimeter Telescope (BLAS T: iDevlin et al.l 120091 ') and 
the Herschel Space Observatory jPilbratt et al. 20ld ) reveal 
strong evolution up to a redshift of at least ~ 0.5, so there 
must be a substantial increase in the numbers and/or lu- 
minosities of dusty star-forming galaxies as we look back in 
cosmic time. In this paper we ask the question: what are 
the properties of dust in galaxies that are not selected to be 
dusty, and is there an evolution in their dust content with 
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redshift equivalent to that seen in Herschel-selected galax- 
ies? 

Galaxies in the Universe comprise an extremely varied 
population, with a wide range of different properties. The 
galaxies that we will concentrate on in this paper are the 
quintessential Hubble tuning fork types, both spirals and el- 
lipticals, that compr ise the majority of galaxies selected in 
optical surveys (e.g. iDriver et alll2006h . We make no prior 
selection with respect to dust content or FIR luminosity, but 
it may be expected that the typical galaxies sampled are qui- 
escent in nature, and are not undergoing excessive starburst 
or nuclear activity (as in typical FIR-selected samples from 
IRAS or Spitzer). This sample may have more in common 
with the low redshift population in the H-ATLAS sample 
selected at 250 /im, which typically consists of optical ly lu- 
minous {M r < -20 ), blue (NUV - r < 4.5) galaxies (pill ; 
iDariush et allhoill ); but unlike H-ATLAS this sample will 



G09 



not be biased towards dusty galaxies in any way. 

Most large statistical studies of the FIR/submm prop- 
erties of FIR-faint galaxies selected by their stellar light 
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Takagi et all 200?l; Serieant et al. 


20081; iGreve et al. 


201C 


Marsden et al. 


120091; Oliver et al. 


2010bl; IViero et al 


I201C 


; iBourne et al 


l201ll). Studies of 



the FIR/submm properties of samples of normal galaxies 
at low/intermediate redshifts have been restricted to small 
sample sizes and most have therefore focussed more on indi- 



Tuffs et al.ll2002l;lLeeuw et al.l2004; Stevens. Amure & 


Gear 
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Stickel, Klaas & Lemke 2007; 


Savov. Welch & Fich 


2009; 


Temi, Brighenti & Mathews! 20091). This is simolv because 



deep submm imaging of large areas of sky is necessary to 
cover a large enough sample of low-redshift galaxies for sta- 
tistical analysis. Until very recently, such data have not been 
available. Observations in the submm, over the Rayleigh- 
Jeans tail of the dust SED at > 200 /im, are crucial for con- 
straining the mass of cold dust in the ISM of galaxies, since 
FIR studies using IRAS at < 100 (jm were only able to con- 
strain the more luminous but less m assive contribution fr om 
warm dust in star-forming regions (|Dunne fc Eales!l200ll ). 

The Hersch el Astrophysical Te rahertz Large Area Sur- 
vey (H-ATLAS; lEales et all l2010bl ) is the first truly large 
area submm sky survey, and as such is ideal for this work. It 
is the largest open-time key project on the Herschel Space 
Observatory and will survey 550 deg 2 in five channels cen- 
tred on 100, 160, 25 0, 350 and 50 ttm, using the P ACS 
(jPoriitsch et all l2010l ) and SPIRE (jGriffin et alj|20ld ) in- 
struments. In this study we use SPIRE maps of the three 
equatorial fields in the Phase 1 Data Release, which cover 
135 deg 2 centred at R.A. of 9 h , 12 h and 14.5 h along the celes- 
tial equator. We are currently unable to use the H-ATLAS 
PACS maps for stacking due to uncertainties in the flux 
calibration at low fluxes, hence this will be pursued in a 
follow-up paper. 

For our galaxy sample we make use of UV, op- 
tical and near-infrared (NIR) data from the Galaxy 
and Mass Assembly (GAMA) survey (|Driver et alj 12009 ) 
which overlaps with the H-ATLAS equatorial fields at 
Deo -1.0° in the 9 h field and Dec> -2.0° in the 
other fields (Fig. [TJ. The GAMA survey combines opti- 
cal data from the Sloan Digital Sky Survey (SDSS DR6; 
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Figure 1. Comparison of the H-ATLAS (blue) and GAMA (red) 
coverage in the three equatorial fields. Dotted lines show the 
H-ATLAS tiles which are not covered in the current data; these 
will be included in a future release. The sample used in this work 
lies within the overlapping area between current H-ATLAS and 
GAMA coverage, approximately 126 deg 2 . 



lAdelman-McCarthv et"alll2008t ), NIR data from the UKIRT 
Infrared Deep Sky Survey (UKIDSS ) Large Area Sur- 
vey (LAS DR4; lLawrence et alj 120071) . and UV from th e 
Galaxy Evolution Explorer {GALEX; Morrissev et al.ll2005l ). 
with redshifts measured with the Anglo-Australian Tele- 
scope and supplem ented with existing redshift surveys (see 
IDriver et al .1120111 for further details) . 

In this paper we select galaxies in the r-band, bin by 
their stellar properties derived from the UV-NIR, and in- 
vestigate their dust properties in the submm. We employ a 
stacking technique to recover median submm fluxes without 
being limited by detection limits in the H-ATLAS images. 
The optical data and sample selection are described in Sec- 
tion [2l and the submm imaging and stacking procedures are 
described in Section although some more technical as- 
pects of the stacking algorithm are left to Appendix [A] In 
Section [4] we present the results of stacking as a function 
of redshift, colour, mass and luminosity, and we discuss im- 
plications for the sources of the extragalactic background. 
In this Section we also derive rest-frame luminosities by in- 
specting the stacked SEDs, and investigate the effects of 
alternative binning schemes on our results. Section [5] con- 
tains a more in-depth discussion of some of the implica- 
tions of the stacking results with respect to dust luminosity, 
temperature and mass. Final conclusions are summarised in 
Section [6] Throughout this work we assume a cosmology of 
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Q A = 0.7, Q M = 0.3, and H a = 70 km s~ x Mpc -1 . All ce- 
lestial coordinates are expressed with respect to the J2000 
epoch. 



2 OPTICAL DATA 
2.1 Sample selection 

We base our selection f unction on the GAMA 'Main Sur- 
vey' (|Baldrv et al.ll2010T ). selecting objects from the GAMA 
catalogue which are classified as galaxies by morphology 
and optical/NIR colours, and are limited in magnitude to 
r-pctro < 19.8 or (z modc i < 18.2 and r modc i < 20.5) or 
(JGnodci < 17.6 and r mo dd < 20.5)0 In fact only 0.3% 
of the sample has r pe tro > 19.8, so the sample is effec- 
tively r selected. To simplify the selection function we use 
the same selection in all fields, so we go below the GAMA 
'Main Survey' cut of r pe tro < 19.4 in the 9 h and 15 h 
fields. For each galaxy we have matched-aperture Kron 
photometry in nine bands: ugriz from SDSS and YJHK 
from UKIDSS-LAS, plus FUV and NUV photometry from 
GALEX More details of the GAMA photometry can be 
found in IHill et all (|201ll l. All magnitudes used in this pa- 
per are corrected for galactic extincti o n usin g the reddening 
data of lSchlegel. Finkbeiner fc Davis! (|l998T l and are quoted 
in the AB system. The Kron magnitudes from IHill et all 
|201ll ) are used for the colours described in this Section, 
however we use the Petrosian measurements for absolute 
magnitudes M r . We purposely do not apply dust correc- 
tions based on the UV-NIR SED or optical spectra, because 
we want to study dust properties as a function of empir- 
ical properties, excluding as much as possible any bias or 
prejudice to the expected dust content. 

We use spectroscopic redshifts from GAMA (year 3 
data) where they are available and reliable (flagged with 
Z_QUALITY (nQ) ^ 3). These are supplemented with pho- 
tometric redshifts computed from the optical-NIR photom - 
etry using ANNZ (for more details see ISmith et alTl2011a[ ). 
The comparison of photometric and spectroscopic redshifts 
is shown in Fig. [5] For this work we apply an upper limit in 
redshift of 0.35, because the number of good spectroscopic 
and photometric redshifts drops rapidly beyond this point. 
This means that the redshift bins would have to be made 
wider to sample the same number of sources (hence dilut- 
ing any sign of evolution); it also means that we only sam- 
ple the brightest objects at the highest redshifts and their 
properties may not be comparable to the typically fainter 
objects sampled at lower redshift. We use photometric red- 
shifts with relative errors Sz/z < 20% only, which excludes 
most of the poor matches in Fig. [2] In this way we exclude 
eight per cent of the whole sample at z < 0.35 (seven per 
cent of Zphot < 0.35), which means that the sample is still al- 
most completely r-band limited. The limiting redshift error 
translates to a 20% error on luminosity distance, a 40% er- 
ror on stellar mass, and an absolute magnitude error of 0.3. 
The criterion tends to exclude lower redshift objects, lead- 
ing to a relative paucity of photometric redshifts at z < 0.2. 
This does not pose a problem since there is near complete 

2 Model magnitudes are the bes t fit of an exponen tial and a 
de Vaucouleurs fit as described bv lBaldrv et al. I i20ld ). 
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Figure 2. Comparison of spectroscopic (nQ ^ 3) and photomet- 
ric redshifts for the objects in our sample which have both. Black 
points have photometric redshifts with relative errors < 20% 
while orange points have greater errors. Using this limiting error 
and a limiting redshift of 0.35 (dashed lines) ensures a reliable set 
of photometric redshifts for our purposes. 

spectroscopic coverage at these lower redshifts. Overall, 90 
per cent of the redshifts we use are spectroscopic, although 
the photometric fraction does increase with redshift out to 
2 = 0.35. The histograms of spectroscopic and photometric 
redshifts are shown in Fig. [3] We tested the effect of random 
photometric redshift errors on our results by perturbing each 
photometric redshift by a random shift drawn from a Gaus- 
sian distribution with width a = Sz. After making these 
perturbations, we made the same cuts to the sample and 
repeated all the analysis, and found that all stacked results 
were robust, changing by no more than the error bars that 
we show. 

We note that a substantial number of photometric red- 
shifts at z > 0.3 appear to be biased low in Fig. [2] This 
explains why there appear to be more photometric redshifts 
than 'all' redshifts at 0.3 < z < 0.35 in Fig. [3]- i.e. some of 
those objects have spectroscopic redshifts which are greater 
than 0.35 and hence do not appear in the same bin in the 
'all redshifts' histogram. This issue could potentially affect 
the results in our highest redshift bin (z > 0.3); the ef- 
fect would be to contaminate that bin with galaxies from 
a slightly higher redshift, which may complicate any evolu- 
tionary trends seen across the z = 0.3 boundary. We have 
chosen to leave the bin in our analysis because over 70 per 
cent of its galaxies have reliable spectroscopic redshifts, so 
the effect of a biased minority of photometric redshifts is 
considered to be small (and ultimately none of our conclu- 
sions hinge on this bin alone). 

In total we have a sample of 86,208 optically se- 
lected galaxies with good spectroscopic or photometric 
redshifts within the 126 deg 2 overlapping area of the 
SPIRE masks and the GAMA survey. We calculated k- 
correc tions for the UV-NIR p hotometry using KCORRECT 
v4.2 (|Blanton fc Roweisl 120071 1. with the spectroscopic and 
photometric redshifts described above. The final compo- 
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Figure 3. Histogram of rcdshifts available in the catalogue: spec- 
troscopic with nQ > 3 (black line); photometric with Sz/z ^ 0.2 
(red line); all rcdshifts combined (grey shaded bars). In construct- 
ing the grey histogram we take all the spectroscopic rcdshifts in 
the black histogram and add any additional photometric redshifts 
from the red. 



nent of the input catalo gue is the set of stellar masses 



fro m iTavlor et all (I201 ll), which were computed by fit- 
ting iBruzual fc Chariot! (|2003l ) stellar pop ulation models to 
the GAMA ugriz photometry, assuming a IChabrierl (|2003T ) 
IMF0 Altogether we have stellar masses estimated for 90 
per cent of the sample. The reason that 10 per cent are 
missing is that our sample reaches deeper than the GAMA 
Main Survey in two of the fields: we use the same magni- 
tude limit in all three fields so that we can sample as large 
a population as possible. 



2.2 Colour classifications and binning 

A simple way to divide the sample in terms of stellar prop- 
erties is to make a cut in rest-frame optical colours. The 
bands which have good signal-to-noise data for the whole 
sample are the three central SDSS bands, hence the most 
reliable and complete optical colours to use are g — r, g — i, 
or r — i. We found very little difference between the distribu- 
tions of colours in any of these three alternatives; each ap- 
pears equally effective at defining the populations of galaxy 
colours. We chose to use g — r since these bands have the 
best signal-to-noise hence greatest depths, and we plot the 
colour-magnitude diagram (CMD) in Fig. [4] 

In this Figure the colour-magnitude space is sampled by 
a two-dimensional histogram in which t he number den sity 
in each bin is weig hted by ]Tl/V ma x jSchmidtl 111)681 ) to 
correct for the incompleteness of a flux-limited sample. To 
achieve this we weighted each galaxy by the ratio of the 
volume of the survey (the comoving volume at z — 0.35) to 
the comoving volume enclosed by the maximal distance at 
which that galaxy could have been detected and included 
in the survey. To measure the latter we use the r pctr o limit 
which is the primary limiting magnitude of the sample (a 



3 The NIR photometry was not used in deriving stellar masses 
due to problems fitting the UKIDSS bands as discussed by 



Figure 4. Two-dimensional histogram (CMD) of rest-frame g — 
r and M r for the full sample. Contours mark the logio of the 
histogram density function, weighted with the 1/V ma x method to 
account for incompleteness as described in the text. Contours are 
smoothed by a Gaussian kernel with width equal to 0.8 of the bin 
width (1/50 of the range in each axis). The heavy black line is the 
red sequence fit given by equation (fTJ . The two dashed lines mark 
the boundaries between red/green and grccn/blue classifications 
respectively, which are given by equation pt . 



negligible proportion of sources that were selected by z or 
K have fainter r pc t ro magnitudes). Our redshift limit of 0.35 
was also considered (so no V ma .x is greater than the comoving 
volume at z = 0.35). 

It is now well established that the optical colours 
of galaxies fall into a bimodal distribution featuring a 
narrow 'red seque nce' and a more di spersed 'blue c loud' 



dTresse et al.lll999l;IStrateva et al]|200ll;lBlanton et al 



Kauffmann et al. 2003: Baldrv et all 12004 iBell et all 



2003 



2004 



Faber et al1l2007l . etcl lBaldrv et all (|2004l ) have shown that 
the optical CMD can be successfully modelled as the sum of 
two Gaussian functions in colour, which evolve with absolute 
magnitude and redshift. In Fig. Owe see a clear bimodality 
in (g — r)rest which can be modelled as a function of abso- 
lute magnitude M r (Petrosian) by splitting the distribution 
into eight bins between M r of —15 and —23, and computing 
the one-dimensional histogram of colours in each bin. These 
histograms were each fitted with the sum of two Gaussian 
functions, shown in Fig. For convenience these functions 
can be thought of as representing two populations, one peak- 
ing on the red sequence and one in the blue cloud, although 
this interpretation has limited physical meaning. A linear 
least-squares fit representing the red-sequence was obtained 
from the means and standard deviations of the red popu- 
lation as a function of absolute magnitude across the eight 
bins. Note that the M r —19.1 bin has effectively no red 
sequence, and has no contribution to the linear fit because 
the standard deviations were used as errors in the fitting. We 
checked for any redshift dependency by splitting the popu- 
lation into three redshift bins as well as magnitude bins, but 
since no variation was found we use the red sequence fit to 
the eight magnitude bins with no redshift binning. This fit 
is shown as a heavy black line in Fig. [3] and is given by 



(g - r) rcBt = 0.724 - 0.026(A/ r + 20). 



(1) 



ITavlor et all j201lT) . 



In order to divide red and blue populations as 
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Figure 5. Histograms of rest-frame g — r colours split into eight 
M r bins between —23 and —15 (weighted with the 1/V ma x 
method). Overlaid are the two Gaussian functions which were si- 
multaneously fitted to each histogram, representing the blue and 
red populations respectively, as well as the sum of the functions. 




0.4 0.6 

(g-r)-0.026(Mr + 20) 

Figure 6. Histogram of rest-frame, slope-corrected g — r colours 
across —23 < M r < —18 (weighted with the 1/V max method). 
Overlaid are the two Gaussian functions fitted to the histogram 
as well as the sum of the functions. The dashed vertical lines show 
the boundaries of the colour bins described in the text. 



confidently as possible we examine the distribution of 
g — r colours across the range —23 < M r < —18 in 
Fig. [6] Here we plot the one-dimensional histogram of 
Cgr = (g — r) IC st — 0.026(A/ r + 20), thus removing the slope 
in the red sequence to emphasise the bimodality. In Fig. [5] 
we exclude M r > — 18 because the red sequence becomes 
negligible at these faint luminosities and the distribution 
becomes dominated by the blue cloud, which hinders our 
two-component fitting (note we do not make any absolute 
magnitude cut when stacking). The distribution in Fig. [6] is 
fitted with the sum of two Gaussians: the red sequence has 
a mean of fj, r — 0.71 and standard deviation of a T = 0.09; 
the blue cloud has fib = 0.43 and at, = 0.14. To make a 
clean sample of red galaxies we make a cut at C gr > 0.67 
(i.e. fi T — 0.5(j r ). This cut was chosen to minimise the con- 
tribution of the 'blue' functional fit, while also including the 
majority (55 per cent) of the area under the red fit. The frac- 
tion of this histogram at C gr > 0.67 that belongs to the blue 
function is seven per cent. It is recognised that the functional 
fits do not necessarily represent two distinct populations of 
galaxies, and this fraction does not imply a contamination of 
the red bin since all galaxies with C gr > 0.67 are empirically 
red. These cuts are largely arbitrary and the main purpose 
they serve is to separate the two modes of the colour dis- 
tribution. Using similar arguments, we make a blue cut at 
C gr < 0.57 (i.e. fib + lo"b) which selects 86 per cent of the 
blue function; the fractional contribution of the red function 
in this bin is four per cent. The intermediate bin by its very 
nature is likely have a heterogenous composition including 
some galaxies close to the red sequence, others that are part 
of the blue cl oud, and some proportion of genuine 'green va l- 
ley' galaxies (|Schiminovich et aljfeoOTl ; iMartin et al.ll2007r i. 
The relative contributions from each of these to the inter- 
mediate ('green') bin may vary as a function of redshift and 
absolute magnitude, which must be kept in mind when draw- 
ing any conclusions. However the red and blue bins will be 
dominated by completely different populations with respect 
to each other at all redshifts and absolute magnitudes, which 
justifies the use of an intermediate bin to fully separate 
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Figure 7. Two-dimensional histogram (CMD) of rest-frame 
NUV — r and M r for the full sample. Contours mark the logio of 
the histogram density function, again using the 1/Vma.x method 
to account for incompleteness as described in the text. Contours 
are smoothed by a Gaussian kernel with width equal to one bin 
(1/40 of the range in each axis). The noise in the top right is 
a result of a small amount of data close to the NUV detection 
limit having exceptionally high weights. The heavy black line is 
the red sequence fit to our data given by equation ([3j . The two 
dashed black lines mark the boundaries between red/green and 
green/blue classifications respectively, which are given by equa- 
tion (3). 

them. The g — r colour bins are summarised in equation ((2j: 

RED: 0.67 + /(M r ) < (g - r) lesi < 1.00 
GREEN: 0.57 + f(M r ) < (g - r) ICBt < 0.67 + f(M r ) 
BLUE: 0.00 < (g - r) rcs t < 0.57 + f(M r ) 

where f(M r ) = -0.026(M r + 20) = -0.026M r - 0.52 (2) 

These divide the sample into 41,350 blue, 17,744 green and 
27,114 red galaxies. 

A limitation of optical colours such as g — r is the small 
separation in colour space between the red and blue popula- 
tions. This would be improved by using a pair of bands which 
straddle the 4000A spectral break, but the only Sloan band 
blueward of this is u, which has poor signal-to-noise and 
therefore a relatively shallow magnitude limi t . It has been 
shown in the literature (e.g . lYi et all 120051 : IWvder et al.l 
l2007l : ICortese fc Hughe J2009r ) that a UV-optical colour such 
as NUV — r (or UV-NIR such as NUV — H) provides greater 
separation between red and blue and reveals a third popula- 
tion of galaxies in the green valley. A clear delineation of the 
populations would minimise contamination between the bins 
and should help to disambiguate trends in stacked results. 

In Fig. [7] we plot the NUV - r CMD (since the r-band 
data are deeper than H , and this gives an iVCTV-limited sam- 
ple). We are unable to successfully fit the NUV — r colour 
distribution in bins of M r using the simple double Gaus- 
sian function, due to a s ignificant 'green va lley' excess. To 
overcome this we follow IWvder et all (|2007h by defining a 
clean red sample [NUV — r > f(M r ) — 0.5] and blue sample 
[NUV - r < /(M r ) - 2.0], where f(M r ) = 1.73 - 0.17M r 
is the fi t to a morphologically-selected red sequence by 
lYi et all (|2005l ). To these subsets we attempt to fit Gaus- 



sian functions for the red and blue distributions respectively, 
again in eight bins of M r between —15 and —23. As before, 
we find a linear least-squares fit to the means of the red 
sequence, given by 

(NUV - r) ICBt = 5.23 - 0.08(M r + 20). (3) 

Th e slope of this fi t is somewhat flatter than the —0.17 found 
bv lYi et alj (|2005t ). but the uncertainty is large due to the 
fact that we have only performed the Gaussian fit in each 
bin to the upper part of the red sequence. As in Fig. [6] we 
could subtract the slope and plot the histogram of the resid- 
ual colour across all M r , but due to the uncertainty on the 
slope this does not give any extra bene fit. We opt instead to 
simply adopt the colour cuts defined bv lWvder et al.l (|2007l ) 
as boundaries between our blue, green and red samples: 

RED: /(M r )-0.5 < (NUV - r) rost < 7.0 
GREEN: /(M r )-2.0 < (NUV - r) rcst < f(M r ) - 0.5 
BLUE: TO < (NUV - r) rcst < f(M r ) - 2.0 

where f(M r ) = 1.73 - 0.17M r (4) 

These divide the sample into 36,900 blue, 12,758 green and 
3115 red galaxies. These numbers reveal two disadvantages 
of using the NUV — r colour: that there are NUV detections 
for only about 60 per cent of the sample, and that the NUV 
selection naturally disfavours red NUV — r colours leading 
to a smaller number of galaxies in the red bin. In contrast 
the r-band selection of the g — r sample is relatively unbiased 
by the colour of the galaxy. However the differentiation be- 
tween blue and green populations should be more successful 
using NUV — r compared with g — r. Therefore both alter- 
natives have their merits. In Section 14.5.11 we compare the 
results obtained using the two alternative colour cuts, but 
for the bulk of the paper we refer to the g — r colour cuts 
unless otherwise stated. 

In this paper we do not explicitly attempt to distinguish 
passive red galaxies from obscured, star- forming red galax- 
ies; rather we focus on the submm properties as a function of 
observed optical colours. We therefore may expect a some- 
what mixed population in the red (and green) bins, even 
using NUV — r. There are various ways one might attempt 
to overcome this - applying dust corrections based on UV 
photometry or spectral line indices, or using the Sersic index 
to predetermine the expected galaxy 'type' - however we 
opt to avoid biasing any of our results by anyprior assump- 
tion about the nature of galaxies in each bin|j We leave any 
analysis that accounts for morphology or spectral properties 
to a follow-up paper. 

3 SUB-MILLIMETRE DATA AND STACKING 

3.1 Stacking into the SPIRE maps 

For the submm imaging we use SPIRE images at 250, 
350 and 500 fim of the three equatorial GAMA fields in 

4 The exception to this empirical approach is the assumption of 
the fc-correction, which is necessary to make fair comparisons be- 
tween redshift bins. The fc-corrections are very well constrained 
by photometry in 5—11 bands and uncertainties are small com- 
pared with our colour bins. 
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H- ATLAS, which were made in a similar way to the s cience 
demonstration maps described bv lPascale etail (|201ll ). The 
fields consist of 53.25 deg 2 at 9 h , 27.37 deg 2 at 12 h 53.93 deg 2 
at 14. 5 h . Background subtraction was c a rried out using the 
NEBULISER routine developed by llrwinl |2010l ') which effec- 
tively filters out the highly varying cirrus emission present 
in the H- ATLAS maps, as well as extended background emis- 
sion including the Sunyaev-ZePdovich effect in clusters and 
unresolved clustered sources at high redshift. 

All sources are treated as point sources and flux den- 
sities (hereafter 'fluxes') are measured in cut-outs of the 
map around each optical position, convolved with a point- 
response function (PRF). We account for sub- pixel scale po- 
sitioning by interpolating the PRF from the point-spread 
function (PSFj] at a grid of pixels offset from the centre by 
the distance between the optical source centre and the near- 
est pixel centre. This convolution with the PRF is the stan- 
dard technique to obt ain the minim um-variance estimate of 
a point source's flux (|Stetsonlll987l '). The PSFs at 250, 350 
and 500 /im have full-widths at half maximum (FWHMs) 
equal to 18, 25 and 35 arcsec respectively. 

We then measure fluxes in the maps at the positions 
in the optical catalogue, and use simplifying assumptions to 
account for blended sources which would otherwise lead to 
over-estimation of flux in the stacks. The procedure is de- 
scribed in detail in AppendixfX] This automatically corrects 
stacked fluxes for the effect of clustered sources in the bins, 
with the caveat that sources not in the catalogue (i.e. be- 
low the optical detection limits) cannot be deblended. We 
stack fluxes by dividing the sample into three colour bins (as 
described in the previous Section), then splitting each into 
five bins of redshift, then six bins of absolute magnitude. 
Redshift and magnitude bins are designed to each contain 
an approximately equal number of objects; in this way we 
ensure that the sample is evenly divided between the bins to 
maintain good number statistics in each. All three fields are 
combined in each stack. We choose to use the median statis- 
tic to represent the typical flux, since the median value with 
a suitable error bar is a robust representation of the distribu- 
tion of fluxes in a given bin. Unlike the mean, the median is 
insensitive to bias due to exceptionally bright s ources which 
are o utliers from the general population (e.g. IWhite et al.l 
l2007h . 

We also measure the background in the maps, since al- 
though they have been sky-subtracted to remove the highly 
variable cirrus foreground emission, the overall background 
does not average to zero, and therefore has a significant con- 
tribution to stacked fluxes. In each map we create a sample 
of 100,000 random positions within the region covered by 
our optical catalogue, masking around each source in the 
main stacking catalogue with a circle of radius equal to the 
beam FWHM in order to avoid including the target posi- 
tions in the background stack. We then perform an identical 
stacking analysis at these sky positions, but reject any po- 
sitions where we measure a flux greater than 5a. This pre- 
vents a bias of our background measurement from sources 



5 The important distinction between the PSF and the PRF is 
that the PRF represents the discrete response function of the 
detector pixels to the continuous distribution of light (PSF) which 
reaches the detector from an astronomical point source. 



detected in the SPIRE maps which have not been masked 
because they are not in the GAMA catalogue. The stacked 
flux measured in this way is a reasonable estimate of the 
average background flux in the corresponding map, and is 
subtracted from our stacked fluxes prior to further analy- 
sis. The average background levels are 3.5 ± 0.1, 3.0 ± 0.1 
and 4.2 ±0.2 mjy beam -1 in the 250, 350 and 500 /xm bands 
respectively. 

Fluxes measured in the SPIRE maps are calibrated for 
a flat vS v spectrum (S„ tx v~ L ), whereas thermal dust emis- 
sion longward of the SED peak will have a slope between t>° 
and v 2 depending on how far along the Rayleigh- Jeans tail a 
given waveband is. The SPIRE Observers' Manua0 provides 
colour corrections suitable for various SED slopes, including 
the v 2 slope appropriate for bands on the Rayleigh- Jeans 
tail. This is a suitable description of the SED in each of 
the SPIRE bands at low redshift, and we therefore mod- 
ify fluxes by the colour corrections for this slope: 0.9417, 
0.9498 and 0.9395 in the three bands respectively. At in- 
creasing redshifts, however, a cold SED can begin to turn 
over in the observed-frame 250 /im band. From inspection of 
single-component SEDs fitted to our stacks we estimate that 
the SPIRE points in most of our bins fall on the Rayleigh- 
Jeans tail, although at the highest redshifts slopes can be 
as flat as v° at 250 /im and v 15 at 350 /im. The correspond- 
ing colour corrections are 0.9888 and 0.9630 respectively. We 
tried applying these corrections to the highest redshift bins 
and found no discernable difference to any of our stacked 
results, hence our results are not dependent on the colour 
correction assumed. 

3.2 Simulations 

The stacking procedure that we used was tested on sim- 
ulated maps to ensure that we could accurately measure 
faint fluxes when stacking in confused maps with realistic 
noise and source density. As described in Appendix [A] wc 
were able to accurately reproduce median fluxes and correct 
errors, although fluxes of individual sources could be under- 
/over-estimated if they were blended with a fainter/brighter 
neighbouring source. 

In addition we simulated various distributions of fluxes 
to test that the median measured flux is unbiased and rep- 
resentative when stacking faint sources close to and below 
the noise and confusion limits. The results of these simula- 
tions indicate that the median can be biased in the presence 
of noise (see also IWhite et all 120071 1. We show details of 
the simulations in Appendix [B] In summary, we assume a 
realistic distribution of fluxes described by dN/dS oc S~ 2 , 
Smin < S < S max ; R = log 10 (S m ax/S'min) = 1-3, and for this 
we estimate the amount of bias in the measured median as 
a function of the true median flux, and correct our stacked 
fluxes for this bias. Correction factors are all in the range 
0.6 — 1.0, and the effect is greatest for low fluxes (< 10 mjy). 
If we consider the true median to be representative of the 
typical source in any bin, then relative to this, the bias in the 
measured median is always less than or equal to the 'bias' 
in the mean resulting from extreme values (as we explain 
in Appendix [Bj. We tested the sensitivity of the results to 

6 Available from http : //herschel . esac . esa. int/Docs/SPIRE/html/spire_om. htm] 
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assumptions about the flux distribution, and found that al- 
though the level of bias does depend on the limits and slope 
of the distribution, all of our measured trends remain sig- 
nificant and all conclusions remain valid for any reasonable 
choice of distribution. This is equally true if we make no 
correction to the measured median. 

We also tested the correlations found in the data by 
simulating flux distributions with various dependencies on 
redshift, absolute magnitude and colour. We first made sim- 
ulations in which fluxes were randomised with no built-in 
dependencies but with the same scatter as in the real data, 
and saw that stacked results were equal in every bin (as ex- 
pected). We tried simulations in which fluxes varied with 
redshift as a non-evolving LF (simply varying as the square 
of the luminosity distance, modified by the fc-correction) , 
with realistic scatter; and also as an evolving LF (log flux 
increasing linearly with redshift at a realistic rate), also with 
scatter. We also allowed flux to vary linearly with optical 
colour and logarithmically with M r , again including realis- 
tic scatter. In all simulations we found that we were able to 
recover the input trends by stacking. 

Finally we simulated fluxes in the three SPIRE bands 
to produce a randomised distribution of submm colours 
(<S25o/£?35o etc) with scatter similar to that in the data but 
with no correlations built in. Results showed that no artifi- 
cial correlations were introduced by stacking. 

These results indicate that the correlations detected in 
the real data (described in the following chapter) should 
be true representations of the intrinsic distributions in the 
galaxy population, and are not artifacts created by the 
stacking procedure. 

3.3 Errors on SPIRE fluxes 

Errors on stacked fluxes are calculated in two different ways. 
Firstly we estimate the instrumental and confusion noise 
on each source and propagate the errors through the stack- 
ing procedure. We estimate the instrumental noise by con- 
volving the variance map at the source position with the 
same PRF used for the flux measureme nt. To this we add i n 
quadrature a confusion noise term, as in lRigbv et al.l (|201 ll ). 
Since fluxes are measured by filtering the map with a ker- 
nel based on the PSF (see Appendix [S} , we need to use 
the confusion noise as measured in the PSF-filtered map. 
iPascale et ail (|201lT ) measured confusion noises of 5.3, 6.4 
and 6.7 mjy per beam in the PSF-filtered H-ATLAS Sci- 
ence Demonstration Phase (SDP) maps. We estimate that 
the confusion noise is at a similar level in our maps after 
PSF-filtering, by comparing the total variance in random 
stacks on the background to the average instrumental noise 
described above. Hence we combine these values of confusion 
noise with the measured instrumental noise of each source. 
In each stack the mean of these measured noises, divided 
by the square root of the number of objects stacked, and 
combined with the error on background subtraction, gives 
the 'measurement error' (on) in Table [1] 

We also use a second method to calculate errors on 
stacked fluxes (as well as other stacked q uantities), which i s 



the 1<7 erro r on the med i an desc ribed bvlGott et al 



as used in iDunne et al.l (|2009al ) and I Bourne et al 



2001), 



2011) 



errors as well as genuine variation within the bin result- 
ing from the underlying population from which it is drawn. 
Briefly, the method sorts the N values in a bin, assigning 
each a unique rank r between and 1. In the limit of large 
N, the expectation value of the true median of the popu- 
lation sampled is (r) = 0.5, and its standard deviation is 
(r 2 — {r) 2 ) 1 / 2 = 1/V4JV. If the measurement at rank r is 
m(r), then the median measurement is m(0.5), which gives 
the expectation value of the true median of the population 
sampled. The error on this expectation value is then 



m(0.5 + l/-\ 



m(0.5 



This is calculated from the distribution of flux values in the 
bin and so automatically takes into account measurement 



(5) 

This formula gives the 'statistical error' (as) in Table [1] 
These values are typically three to four times larger than 
the measurement error on, indicating that the uncertainty 
resulting from the spread of intrinsic fluxes in a bin is greater 
than the combined noise in the map at all the positions in 
the stack. 



4 RESULTS 

4.1 Stacked fluxes 

The results of stacking SPIRE fluxes as a function of g — r 
colour, redshift and absolute magnitude M r are given in 
Table [T] Secure detections are obtained at 250 and 350 fim 
in most bins, although many bins have low signal-to-noise 
at 500 /im. Note that the signal-to-noise ratios in the Table 
are based on the measurement error reduced by y/N (i.e. 
on), since this strictly represents the noise level (instrumen- 
tal plus confusion) which we compare our detections against. 
When talking about errors in all subsequent analysis we will 
refer to the statistical uncertainty on the median (<rs) be- 
cause this takes into account both instrumental noise and 
the distribution of values in the bin, both of which are con- 
tributions towards the uncertainty on the median stacked 
result. 

Table [T] also contains the results of Kolmogorov- 
Smirnov (KS) tests which were carried out on each stack 
to test the certainty that the stacked flux represents a sig- 
nal from a sample of real sources and is not simply due to 
noise. This was done by comparing the distribution of mea- 
sured fluxes in each bin with the distribution of fluxes in the 
background sample for the same SPIRE band. These back- 
ground samples (described in Section f3.1|l were placed at a 
set of random positions in the map, after masking around 
the positions of input sources. If any of our stacks did not 
contain a significant signal from real sources then the KS 
test would return a high probability that the distribution 
of fluxes is drawn from the same population as the back- 
ground sample. The great majority of our bins were found 
to have an extremely small KS probability, meaning that 
we can be confident that the signals measured are real. The 
highest probability is 0.04, for a 500 fim stack in the highest- 
redshift, red-colour bin. A small sample of the bins are ex- 
plored in more detail in Appendix[Cl where we show stacked 
postage-stamp images and histograms on which the KS test 
was carried out. 

Stacked fluxes at 250, 350 and 500 /im in the observed 
frame are plotted in Fig. |8j showing the dependence on M r 
and g—r at different redshifts. The majority of the bins have 
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stacked fluxes well below the 5a point-source detection lim- 
its shown as horizontal lines on the Figure. In all three bands 
there is a striking difference between the submm fluxes of 
blue, green and red galaxies, and a strong correlation with 
M r in the low-redshift bins of blue and green galaxies. These 
trends unsurprisingly indicate that the red galaxies tend to 
be passive and have lower dust luminosities than blue, and 
are generally well below the detection threshold in all SPIRE 
bands. They also show that submm flux varies little across 
the range of M r in red galaxies, while in blue galaxies it cor- 
relates strongly with M r such that only the most luminous 
r-band sources ha ve fluxes above the 25 /im detection limit 
- a point noted bv lDariush et all (|201lh and ID 111 The vari- 
ation with redshift is also very different between the three 
optical classes, with the fluxes of blue galaxies diminishing 
with redshift more rapidly than those of red galaxies. Fluxes 
in the green bin initially fall more rapidly with increasing 
redshift than those in the blue, but at z > 0.18 they resemble 
those in the red bin and evolve very little. This is potentially 
due to a change in the nature of the galaxies classified as 
green at different redshifts, which is unsurprising since this 
bin contains a mixture of different galaxy types in the over- 
lapping region between the blue cloud and red sequence. It 
is likely that the relative fractions of passive, relatively dust- 
free systems and dusty star-forming systems in this bin will 
change with redshift, as the star-formati on density of the 
Universe evolves over this redshift range l|Lillv et al. 1 1 19961 ; 
iMadau et al.|[l99r3 . etc). The evolutionary trends discussed 
in this Section can be better explored by deriving submm 
luminosities, which first requires a model for fc-correcting 
the fluxes, as we will discuss in Section 14.41 

At this point it is worth considering some potential 
sources of bias in different bins in case they might impact on 
the apparent trends. For example, it is reasonable to expect 
that certain classes of galaxy are more likely than others to 
inhabit dense environments: in particular redder galaxies, 
and more massive galax ies, are known to be more clustered 
(e.g. IZehavi et aT]|201lh . While we do account for blending 
in the flux measurements (see Appendix \K§ this is limited 
to the blending of sources within the catalogue. As we move 
to higher redshifts the catalogue becomes more incomplete, 
and it becomes more likely that the clustered galaxies will 
be blended with some unseen neighbour. We make no cor- 
rection for this effect, but we expect the contamination to be 
small for the following reasons: The input sample is complete 
down to below the knee in the optical LF at z < 0.3 (see next 
section); we therefore account for the blending with most of 
the galaxies in the same redshift range. Contaminating flux 
would have to come from relatively small galaxies which are 
not likely to contribute a large amount of flux. Moreover, 
the blending corrections are on average very small in com- 
parison with the trends that we observe (see Appendix |A1 
Table |AT]) . so a small additional blending correction should 
not significantly alter our conclusions. 

4.2 Contamination from lensing 

There is a risk that the submm fluxes of some galaxies in 
the sample are contaminated by flux from lensed background 
sources, via galaxy-galaxy lensing. This is especially likely in 
a submm survey as a resul t of the negative fc-correction an d 
steep evolution of the LF (|Blainll 19961 : iNegrello et al.ll2007h . 



These factors conspire to make submm sources detectable 
up to very high redshifts, therefore providing an increased 
probability for one or more foreground galaxies to intrude in 
the line of sight and magnify the flux via strong gravitational 
lensing. The strong potential for detecting lensed high red- 
shift sources in H-AT LAS was conclusively demonstrated by 
INegrello et all (|2010h . In this study we target low redshift 
sources selected in the optical, but our sample will inevitably 
include some of the foreground lenses whose apparent fluxes 
are likely to be boosted by flux from the lensed background 
sources. The flux magnification is likely to be greatest for 
massive spheroidal lenses, a s a result of their mass distri- 
bution (|Negrello et alj|20ld ). This could pose a problem for 
our red bins, where spheroids will be mostly concentrated. 
To make matters worse, measured fluxes are lowest in our 
red bins which means that even a small lensing contamina- 
tion of order 1 mjy could significantly boost the flux. 

We can make an estimate of the lensing contribution to 
stacked fluxe s by considering t he predicted number counts 
of lenses from lLapI et al.l(|201ll 'l. which are based on the am- 
plification distribution of s t rong l enses (amplification factors 
> 2) from INegrello et"aH (|2007T ). Integrating these counts 
gives a total of 470 lensed submm sources per square de- 
gree, and integrating their fluxes per square degree gives 
the total surface brightness of lensed sources shown in the 
first line of Table [2] However, the counts are not broken 
down by redshift, and only those at z < 0.35 will contribute 
to our stacks. It is not trivial to predict what fraction of 
strong lenses are in this redshift range, but recent results 
from H-ATLAS can provide us with the bes t estim ate that 
is currently possible. iGonzalez-Nuevo et al.l |2012f ) created 
a sample of 64 candidate strong lenses from the H-ATLAS 
SDP by selecting sources with red SPIRE colours which have 
no reliable SDSS IDs, or have SDSS IDs with redshifts incon- 
sistent with the submm SED. After matching to NIR sources 
in the VISTA Kilo-degree INfrared Galaxy survey (VIKING; 
Sutherland, W. et al. in preparation), they reduced this sam- 
ple to 33 candidates with photometric redshifts for both the 
lens (using the NIR photometry) and source (using SPIRE 
and PACS photometry). This sample, the H-ATLAS Lensed 
Objects Selection (HALOS), is unique in being selected in 
the submm, enabling the selection of candidate lenses over a 
much larger redshift range than other lens samples to date 
(their lenses had photometric redshifts < 1.8, while other 
surveys were confined to z < 1). HALOS therefore provides 
the best observational measurement of the lens redshift dis- 
tribution. 

Sev en of the 33 HALOS ca ndidates have lens redshifts 
< 0.35. IGonzalez-Nuevo et al.1 removed two of these from 
their final sample because the lenses were at z < 0.2, and 
they considered lenses at such low redshifts to have low 
probability both on t heoretical grounds and on the evidence 
of previous surve ys (|Browne et al.l [20031 : lOguri et al.l 120061 : 
iFaure et al.ll200St ). However, we do not want to risk under- 
estimating the number of low redshift lenses, so we conser- 
vatively include those two in our analysis. The fraction of 
lenses at z < 0.35 is therefore 7/33 = 21±i? %, using bi- 
nomial techniqu es to estimate the 1-a confidence interval 
(|Cameronll201ll ). Using these results to scale the total lensed 
flux from all redshifts, we obtain the contribution from lenses 
at z < 0.35, as shown in Table [2] Assuming that all these 
low redshift lenses fall in the red bin of our sample, we can 
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Table 1: Results of stacking in bins of g — r colour, redshift, and absolute magnitude (M r ). Columns are as follows: colour bin C — B (blue), G (green), R (red); median 
redshift z in bin (approximate z bin boundaries are 0.01, 0.12, 0.17, 0.22, 0.28, 0.35); median M r \ count N in the bin. The following columns for each of the three SPIRE 
bands: background-subtracted flux S in mjy; signal-to-noise ratio S/a^; measurement error <7n (mjy) computed f rom the mean vari ance of positions in the stack, divided 
by \/~N, and including the error on background subtraction; statistical error on the median flux (as, mjy) following lGott et al.l (|200l[ ); KS probability that the distribution 
of fluxes in each bin is the same as that at a set of random positions; median fc-correction K' — K(z)/(1 + z) in the bin. 
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Figure 8. Stacked SPIRE fluxes (not k— corrected) as a function of g — r colour, redshift and absolute magnitude M r . Top: 250 (im, 
middle: 350 fim, bottom: 500 fim. Galaxies are binned by optical colour from blue to red (shown in panels from left to right) and by 
redshift (shown by plot symbols) and stacked fluxes in each bin are plotted against M r . Error bars are the statistical la errors in the 
bins as described in Section |3 .31 and also include errors due to background subtraction. The horizontal dashed lines at 33.5, 39.5 and 
44.0 mjy in 250, 350 and 500 fim, represent the 5<r point-source detection limits as measured in the PSF-convolved Phase 1 maps. 



compare these fluxes to the total stacked flux of our red bins 
as shown in Table [2] which indicates that about 10, 20 and 
30 per cent of the 250, 350, and 500 fim fluxes respectively 
comes from high redshift sources lensed by the targets. This 
may be a slight overestimate si nce some of t he len ses may 
fall in the other bins; however lAuger et al.l (|2009l ) showed 
that 90 per cent of lenses are massive early type galaxies. 
Any lensing contribution to the blue or green bins would be 
negligible compared to the fluxes measured in those bins. 



The lensed flux is divided between the redshift bins of 
the red sample in a way that is determined by the prod- 
uct of the lens number distribution rii(z) and the lens ef- 
ficiency distribution $ ( z). T he numbers m(z) are given by 
iGonzalez-Nuevo et al.l (|2012l ). while the efficiency depends 
on the geometry between source, lens and observer. We esti- 
mate $(z) from the HALOS so urce and le ns redshift dis- 
tributions using the formula of [Hu| (|l999h . and compute 
the lens flux distribution from the product of total lensed 
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Tabl e 2. Total surface brightness of lensed sources from the 
lLapi et al.l teOllT) counts model, and estimated contribution from 
the low redshift populatio n of lenses assuming the lens redshift 
distribution from HALOS dGonzalez-Nuevo et ai1l2012l) . This is 
compared to the total surface brightness of red galaxies (g — r 
colour) at z < 0.35 from our stacks. We then estimate the frac- 
tion of the flux in each redshift bin of the red sample that comes 
from lensed background sources. 





250 Atm 


350 /im 


500 fan 


All lensed flux 


Total surface brightness 
1.09 1.34 


Jy deg" 2 
1.22 


Lenses at z < 0.35 
Red galaxies 


n 9 o+0.09 
u ' zo -0.06 

2.6 ±0.5 


0.28i°- 
1.6 ±0.2 


26+ 11 
0.8 ± 0.1 


0.01 < z < 0.12 


Lensed flux/red galaxy flux by z bin 
0.00 ±0.00 0.00 ±0.00 0.00 ±0.00 


0.12 < z < 0.17 


0.06 ±0.01 


0.12 ±0.02 


0.21 ±0.04 


0.17 < z < 0.22 


0.11 ± 0.02 


0.20 ±0.04 


0.33 ±0.06 


0.22 < z < 0.28 


0.16 ± 0.03 


0.27 ±0.05 


0.44 ± 0.08 


0.28 < z < 0.35 


0.20 ±0.04 


0.35 ±0.07 


0.58 ±0.11 



flux, ni(z) and $(z). Comparing this to the total flux of red 
galaxies in each redshift bin, we compute the fractional con- 
tamination from lensed flux as shown in Table [2] Errors on 
the lensed flux per redshift bin are dominated by the Pois- 
son error on the normalisation of n;(z), which is simply the 
Poisson error on the count of 33 lenses. The relative error on 
the lensed flux is therefore \/33/33 = 0.17. The error on the 
stacked red galaxy fluxes is dominated by the 7% flux cal- 
ibration error (Pascale E., et al. in preparation), hence the 
errors on the fractions in Table [2] are given by the quadra- 
ture sum of 7% and 17%, which is 19%. Using the fractions 
derived above we can remove the estimated lensed contribu- 
tion to stacked fluxes in each redshift bin of the red sample. 
The effect of subtracting this fraction from the fluxes of red 
galaxies is minor in comparison to the trends described in 
Section f4.ll The effect on other derived results will be dis- 
cussed later in the paper. 

4.3 Resolving the cosmic IR background 

A useful outcome of stacking on a well-defined population 
of galaxies such as the GAMA sample is that we can eas- 
ily measure the integrated flux from this population and 
infer how m uch it contributes to the cosmic infra red back- 
ground (CIB; |Puget et al.ll99tj ; lFixsen et al.lll998T ). The cos- 
mic background at FIR/submm wavelengths makes up a 
substanti al fraction of the integrated radiative energy in the 
Universe (|Dole et al.ll2006t ). although the sourc es of this ra- 
diation are not fully accounted for. For example. lOliver et al.l 
(|2010al ) calculated that the HerMES survey resolved only 
15 ± 4% of the CIB into sources detected with SPIRE at 
250 /im, down to a flux limit of 19 mjy. A greater fraction 
can be accounted for using P(D) fluctuation analysis to 
reach below the d etection limit of the map: in HerMES, 
iGlenn et al.l |2010r i resolved 64 ± 16% of the 250 /mi CIB 
into SPIRE sources with S250 > 2 mjy. Stacking on 24 /im 
sources has also proved successful, utilising the greater depth 
of 24 /im maps from Spitzer-MIPS to determine source cat- 
alogues for stacking at longer w avelengths. Stacking into 
BLAST, iBethermin et all (|201Ch resolved 48 ± 27% of the 
250 /im CIB into 24 /xm sources with S250 > 6.2 mjy while 



iMarsden et~aH |2009T ) resolved 83 ± 21% into sources with 
S24 > 15 /iJy. However, these BLAST measurements in- 
cluded no corrections for clustering; the authors claimed 
that the effect was negligible, although this observation may 
appear to conflict w i th similar analyses in the literature 
(iNegrello et all 120051 ; ISerieant et allhoOSl ; ISerieantJ |2010| ; 
iBourne et al.ll201ll ). 

Similarly we can stack the GAMA sample to estimate 
what fraction of the CIB at 250, 350 and 500 /im is pro- 
duced by optically detected galaxies at low redshifts. To 
do this we measure the sum of measured fluxes in each bin 
and scale by a completeness correction to obtain the to- 
tal flux of all r < 19.8 galaxies at 2 < 0.35. The correc- 
tion accounts for two levels of incompleteness. The first is 
the completeness o f the original magnitude-limited sample: 
iBaldrv et al.1 (|2010T l estimate that the GAMA galaxy sample 
(after star-galaxy separation) is > 99.9% complete. The sec- 
ond completeness is the fraction of the catalogue for which 
we have good spectroscopic or photometric redshifts (i.e. 
spectroscopic Z_QUALITY ^ 3 or photometric Sz/z < 0.2; 
see Section l2.ip . This fraction is 91.9%; however we have 
only included galaxies with redshifts less than 0.35, which 
comprise 86.8% of the good redshifts. We cannot be sure 
of the redshift completeness at 2 < 0.35 (accounting for 
both spectroscopic and photometric redshifts) so we sim- 
ply assume that we have accounted for 91.9% of these, to 
match the redshift completeness of the full sample^ Finally 
we scale by the fraction of galaxies at 2 < 0.35 that are 
within the overlap region between the SPIRE mask and the 
GAMA survey, which is 72.4%. The combined correction fac- 
tor is 77 = 1/(0.999 x 0.919 x 0.724) = 1.504. The corrected 
flux is converted into a radiative intensity (nW mT 2 sr _1 ) by 
dividing by the GAMA survey area (0.0439 sr). We compare 
this to the CIB le vels expected in the three SPIRE bands 
|Glenn et alJl201Ch - thes e are c alculated by integrating the 
CIB fit from lFixsen et al.1 (| 19981 ) over the SPIRE bands. We 
find that the optical galaxies sampled by GAMA account 
for < 5% of the background in the three bands (see Ta- 
ble [3}. In the Table we also show the percentage of the CIB 
produced by galaxies at z < 0.28, since in this range the cat- 
alogue is complete down to below the knee of the optical LF 
at M * — —21.4 (Petrosian magnitude, h = 0.7; iHill et all 

Soul). 



4.4 The submm SED and k— corrections 

Monochromatic luminosities (WHz" 1 ) at rest-frame 250, 
350 and 500 /im can be calculated using equation ([6]), in 
which S v is the SPIRE flux in Jy, K(z) is the fc-correction 
at redshift 2, and Dl the corresponding luminosity distance 
(m). The (1 + z) on the denominator is the bandwidth cor- 
rection, which together with the fc-correction converts an 
observed 250 /im flux to the flux at rest-frame 250 /im. 

_ 26 4 itDl S v K{z) 



L v = 10" 



(6) 



(l + z) 

Tf-corrections are obtained by assuming that the SED emit- 
ted by dust at a temperature T dust is governed by a greybody 

7 This may be a slight underestimate of the redshift complete- 
ness at z < 0.35, in which case we would overestimate the total 
corrected flux by a maximum of 8.7%. 



© 0000 RAS, MNRAS 000, 000-000 



H-ATLAS/GAMA: dust in optically selected galaxies 15 



Table 3. Total intensities of r pe tr < 19-8 galaxies from stacking 
at 250, 350 and 500 (im, i n com parison to the corresponding CIB 
levels from lFixsen et al | il998h . We show the intensity as a per- 
centage of the CIB for the full stack, and for th e z < 0,28 subse t 
which is complete in M r down to M* = — 21. 4 jffill et al.ll201 if) . 
We also show the contributions of the individual redshift bins and 
g — r colour bins. All contributions from red galaxies have been 
corrected for the lensed flux contamination using the fractions in 
Tabic [2] All errors include our statistical error bars from stack- 
ing, the error on the lensing correction (where applicable) and a 
7% flux calibration error (Pascale, E. et al. in preparation). 



fe.g.lDunne fc Ealesll200ll; I James et al.ll2002l; IPopescu et all 
' 20021; iBlain, Barnard fc Chapman! 120031; iLeeuw et allbOolT 
Hill et al.ll2006l : IParadis. Bernard fc Menvll2009h . For com- 





250 fim 


350 fim 


500 nm 




Intensity nW m 2 sr 


-l 


CIB 


10.2 ± 2.3 


5.6 ±1.6 


2.3 ±0.6 


Total Stack 


0.508 ±0.036 


0.208 ±0.015 


0.064 ± 0.005 


0.01 < z < 0.28 


0.428 ± 0.030 


0.173 ±0.012 


0.054 ± 0.004 






% of CIB 




Total Stack 


4.98 ±0.39 


3.71 ±0.30 


2.79 ± 0.22 


0.01 < z < 0.28 


4.19 ±0.34 


3.08 ±0.26 


2.33 ± 0.19 


0.01 < z < 0.12 


1.57 ±0.17 


1.11 ±0.13 


0.81 ± 0.09 


0.12 < z < 0.17 


0.97 ±0.10 


0.71 ±0.08 


0.53 ± 0.06 


0.17 < z < 0.22 


0.85 ±0.08 


0.63 ±0.07 


0.49 ± 0.05 


0.22 < z < 0.28 


0.80 ±0.08 


0.63 ±0.07 


0.50 ±0.05 


0.28 < z < 0.35 


0.78 ± 0.08 


0.63 ±0.07 


0.46 ± 0.05 


Blue 


3.02 ±0.26 


2.34 ±0.21 


1.67 ±0.15 


Green 


1.12 ±0.10 


0.84 ±0.08 


0.64 ±0.06 


Red 


0.83 ±0.08 


0.64 ±0.07 


0.48 ± 0.05 



of the form B(y, Td ua t) (where B is the Planck function) 
The fc-correction for this SED is given by 

3+/3 „h„ /feT dust _ , 



a hv a /kT du 



(7) 



where v Q is the observed frequency in the 250, 350 or 500 
band, i/ e = (1 + z) v a is the rest-frame (emitted) frequency, 
k is the Boltzmann constant and h the Planck constant. 

Since we cannot fit an SED to individual galaxies we 
instead examine the ratios between stacked fluxes in each 
bin. Fluxes in the red bins are first corrected to remove 
the contribution from lensing as discussed in Section 14.21 
The colour-colour diagrams in Fig. [9] show the resulting 
flux ratios in the observed frame in each of the five red- 
shift bins, alongside a selection of models, which are plot- 
ted with a range of temperatures increasing from left to 
right. We try both a single greybody and a two compo- 
nent model, but there is little to choose between them in 
these colours, since the SPIRE bands are at long wave- 
lengths at which the SED is dominated by the cold dust, 
with little contribution from transiently heated small grains 
or hot dust. We therefore adopt a single component for 
simplicity. The scatter in the data is large, as are the 
errors on the 350/500 )im flux ratio. Moreover, with all 
our data points on the longward side of the SED peak 
we are unable to resolve the degeneracy between the dust 
temperature and the emissivity index (/3). This is shown 
by the close proximity of the (3=1 and f3 = 2 models 
on Fig. [9] which overlap in different temperature regimes. 
With these limitations we are forced to assume a constant 
value of /3 across all our bins. We choose a value of 2.0, 
which has been shown to be realistic in this frequency range 



parison the Planck Collaboration found an average value of 
P = 1.8 ± 0.1 by fitting SEDs to data at 12 Ami-21 cm from 
across the Milky Way ([Planck Collaboration et all l2011al . 
also references therein) |f| 

Under the assumption of a constant f3 (whatever its 
value) the dust temperatures implied by Fig. [9] take a wide 
range of values across the various bins (between 11 and 
22 K for f3 = 2). This is not just random scatter; red galax- 
ies tend towards colder temperatures than blue and green, 
while blue galaxies in some redshift bins show a trend to- 
wards lower temperatures at brighter M r . For the purposes 
of fc-corrections we can estimate the temperature more ac- 
curately by fitting greybody SED models to the three data 
points at the emitted frequencies given by the observed fre- 
quency scaled by 1 + z, using the median redshift in the bin0 
In general one must be careful when using stacked fluxes in 
this way to examine the SED, since when stacking many 
galaxies with different SEDs, the ratios between the stacked 
fluxes can be unpredictable and not representative of the in- 
dividual galaxy SEDs. In this case however we believe we can 
be fairly confident of the results because we bin the galax- 
ies in such a way that we should expect the SEDs within 
each bin to be similar, and so inferred dust temperatures 
and other derived parameters should be accurate. 

The best-fit temperatures range between 12 — 28 K, with 
a median value of 18.5 K. Using j3 = 1.5 instead, the tem- 
peratures are increased by a factor 1.2 — 1.6, ranging from 
13 — 46 K with a median of 23.0 K. We can compare this me- 
dian value to temperatures deriv ed from single-co mponent 
fits in the literature. For example. iDve et all (feoiot ) derived 
a median isothermal temperature of 26 K (/3 = 1.5) for the 
detected population in the H- ATLAS science dem onstration 
data, in agreement with the BLAST sample of IDve et ail 
i|2009h . The value of 26 K is within the range of our tem- 
peratures using /3 = 1.5, and only slightly hi gher than the 
media n. Higher temperatures were found by Hwang et all 
(|2010f i in their PEP/HerMES/SDSS sample of 190 local 
galaxies : they reported median temperatures rising as func- 
tion of IR luminosity, from around 26 K at 10 9 Lq to 32 K 
at 10 11 Lq and 40 K at 10 12 L (/3 = 1.5). These temper- 
atures may be higher because iHwang et all required a de- 
tection shortward of the SED peak (i.e. in an AKARI-FIS 
or IRAS band) for galaxies to be included in their sam- 
ple. Fitting a single greybody to an SED which contains 
both a cold (< 20 K ) and a warm (> 30 K) component 
l|Dunne fc E alcs 200j]) may give results that are not com- 
parable to our s, which fit o nly th e cold component. On 
the other hand. lSmith et all (|2011bl ) fitted greybodies with 
ft = 1.5 to the H- ATLAS 250 /im-selected sample of low red- 



8 A further issue wit h fitting SEDs is tha t /3 may vary with fre- 
quency. For example Paradis ct al. (2009) analysed data on the 
Milky Way from 100 fim to 3.2 mm and showed that (3 was gen- 
erally steeper at 100 — 240 fim than 550 — 2100 /jm. This is an 
effect that we cannot take any account of without many more 
photometric points on the SED, but it could have some effect on 
our fitted temperatures and therefore luminosities. 

9 The temperature fits and trends mentioned here are discussed 
further in Scction l5.ll 
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shift galaxies matched to SDSS, and found a median tem- 
perature of 2 2.5 ± 5.5 K (similar to our result), and unlike 
iHwang et al.l they found no evidence for a correlation with 
luminosity. 

The IPlanck Collaboration et all (|2011bh compiled a 
sample of around 1700 local galaxies by matching the Planck 
Early Release Compact Source Catalogue and the Imperial 
IRAS Faint Source Redshift Catalogue, and fitted SEDs to 
data between 60 — 850 fim using both single-component fits 
with variable /3, and dual-component fits with fixed ,8 = 2. 
In their single-component fits they found a wide range of 
temperatures (15 — 50 K) with median T = 26.3 K and me- 
dian p — 1.2. This median temperature is consistent with 
the Herschel and BLAST results, and the low value of /? is 
likely to be due to the inclusion of shorter wavelength data. 
The authors state that the two-component fit is statistically 
favoured in most cases; these fits indicate cold dust tem- 
peratures mostly between ~ 10 — 22 K, consistent with the 
range in our data. 

In any case we do not necessarily expect to find the 
same dust temperatures in an optically selected sample as in 
a submm selected sample. For the purposes of fc-corrections 
this is relatively unimportant, at least at the low redshifts 
covered in this work. The choice between cold T/high f3 and 
hot T/low /? makes very little difference to monochromatic 
luminosities, as crucially they both fit the data. Likewise 
the range of temperatures has little effect on fc-corrections: 
using the median fitted temperature of 18.5 K in all bins 
gives essentially the same results as the using the temper- 
ature fitted to each bin separately. To remove the effect of 
the variation between models, we carry out all analysis of 
monochromatic luminosities using the median temperature 
of Tdust = 18.5 K and f3 = 2.0 to derive fc-corrections using 
equation J7|| (except where stated otherwise). The implica- 
tions of the fitted SEDs on the physical properties of galax- 
ies in the sample will be discussed in Section [S] First we 
will concentrate on the observational results of the stacking 
which are not dependent on the model used to interpret the 
submm fluxes. 

4.5 Luminosity evolution 

To calculate stacked luminosities we apply equation ((6]) to 
each measured flux and stack the results. The error on the 
stacked value is again calculated using the lGott et"aT1 (|200lh 
method. Note that this method is not the same as applying 
equation ([6]) to the median flux and median redshift of each 
bin, since luminosity is a bivariate function of both flux and 
redshift. Fluxes of sources in the red bin are corrected for 
the fractional contributions from lensing given in Table [2l 
as explained in Section f4. 21 The l-cr errors on these correc- 
tions are included in the luminosity errors. Results in Fig. 1 101 
show a generally strong correlation between luminosities in 
the r-band and all three submm bands, as is the expected 
trend across such a broad range, but the dependence is not 
on M r alone. This becomes obvious when comparing the 
data points with the grey line, which shows the linear least- 
squares fit to the results from the lowest-redshift blue galax- 
ies (the line is the same in each panel from left to right). In 
the blue galaxies, there may be a slight flattening of the 
correlation for the brightest galaxies and/or the higher red- 
shifts, but this effect is much stronger in the green galaxies, 



which are intermediate between the blue and red samples. 
For the red galaxies the correlation disappears entirely but 
for the faintest bin at low redshift. The luminosities of the 
red galaxies all lie below the grey line, showing that red 
galaxies emit less in the submm than blue or green galaxies 
of the same M r , strongly suggesting that they are domi- 
nated by a more passive population than green and blue 
galaxies. These trends are greater than the uncertainties on 
the lensing correction. 

Apart from this colour dependence there is also a signif- 
icant increase in submm luminosity with redshift for green 
and red galaxies of the same r-band luminosity. This evolu- 
tion appears to occur at all M r , without being particularly 
stronger for either bright or faint galaxies, but it is especially 
strong for red galaxies. This may indicate a transition in the 
make-up of the red population, with obscured star-forming 
galaxies gradually becoming more dominant over the pas- 
sive population as redshift increases. Such a scenario might 
be expected as we look back to earlier times towards the 
peak of the universal star-formation history. One problem 
with this explanation is that we might expect an increase in 
obscured star-formation to be accompanied by an increase 
in the dust temperatures at higher redshifts, which we found 
no evidence for in the SPIRE colours (Section 14. 4|) . 

Meanwhile the green sample shows similarities with 
the blue at low redshift and low r-band luminosity, but at 
high redshifts and stellar masses the luminosity dependence 
on M r is flatter and more similar to that of red galaxies. 
This could be due to a shift in the dominant population 
of the green bin, between blue-cloud-like galaxies and red- 
sequence-like galaxies at different redshifts and M r . 

4-5.1 UV-optical versus optical Colours 

Splitting the sample by the NUV — r colour index provides 
a slightly different sampling regime and reduces contami- 
nation between the colour bins because the red and blue 
populations are better separated (see Section f2 . 2 [) . It there- 
fore offers a useful test of the robustness of the results of 
stacking by g — r. Figure [11] shows that stacked 250 /im lu- 
minosities follow the same trends with colour, redshift and 
M r as in the g — r stacks (350 and 500 /im results are simi- 
lar). This supports the interpretation that the three colour 
bins sample intrinsically different populations in terms of 
the dust properties. The red sample in either NUV — r or 
g — r appears to be dominated by passive galaxies at low 
redshifts at least, but the emission from dust increases by a 
factor of around 10 over the redshift range. 

Errors are slightly larger in this sample, particularly in 
the red bin, because we are limited to the 52,773 galaxies 
with NUV detections. Since results appear to be indepen- 
dent of the colour index used, we opt to use the more com- 
plete r-limited sample of 86,208 sources with g — r colours 
for all subsequent analysis. 

4-5.2 Stellar mass versus absolute magnitude 

An alternative to dividing the sample by M r is to use the 
stellar masses w hich were calculated from the GA MA ugriz 
photo metry by iTavlor et all (|201ll ) assuming a IChabrieil 
|2003l ) IMF. Stellar mass is a simple physical property of 
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Figure 9. Colour-colour diagram of the observed-frame SPIRE fluxes. The plot is divided into five redshift bins, in which the data from 
that bin are plotted along with the colours expected from various models as they would be observed at the median redshift in the bin. 
Data are divided into the three g — r colour bins, denoted by symbols and colours, and six M r bins denoted by the size of the data point 
(larger=brighter). Data points and error bars in the red bins include the lensing correction and its uncertainty. Three families of models 
are shown: two consist of isothermal SEDs with either = 1 or /3 = 2, and various dust temperatures; the third is a two-component 
SED with /3 = 2, warm dust temperature T w = 30 K, with a cold/hot dust ratio of 100. Each model is given a range of (cold) dust 
temperatures; the dots along the lines indicate IK increments from 10 K (lower left) to 30 K. Choosing a single component model with 
f) = 2 leads to a range of temperatures between 13 and 22 K. 



the galaxy so may reveal more about intrinsic dependencies; 
on the other hand it depends much more on the models used 
to fit the optical SED than M r , which is only subject to a 
small fc-correction and the assumed cosmology (for a given 
redshift). Relative errors on stellar masses are dominated by 
syste matics, but are small (AlogM at ar ~ 0.1; iTavlor et al] 
2011). We confirmed that our results are robust to these er- 
rors by repeating all analysis after making random pertur- 
bations to the stellar masses, where the size of each pertur- 
bation was drawn from a Gaussian distribution with width 
a = AlogM Btar . No results were systematically affected by 
these perturbations, and random deviations in stacked val- 
ues were smaller than the error bars. 

Figure [12] shows that the results of stacking by g — r 
colour and stellar mass differ slightly from the results of 
stacking by M r (Fig. I lOf) . Again we found very little differ- 
ence from these results when we stacked by NUV — r colour 



and stellar mass. The results of stacking by mass seem to 
differ most in the blue bin. Whereas there was little lumi- 
nosity evolution at fixed M r , these results show evolution 
at fixed stellar mass. Furthermore this evolution is depen- 
dent on stellar mass, suggesting that smaller galaxies tend 
to evolve more rapidly. The samples in Figures [10] and [12] 
are slightly different since stellar masses were only avail- 
able for 90 per cent of the full sample; however we know 
that this is not responsible for the discrepancy since repeat- 
ing the stacking by M r with the stellar mass sample gives 
identical results (the stellar mass incompleteness does not 
vary between bins). The difference arises because M r does 
not directly trace stellar mass, which leads to a mixing of 
galaxies of different masses within a given M r bin. This is 
unavoidable since we must split the sample in three ways 
(colour, redshift and mass/magnitude) because dust lumi- 
nosity varies strongly as a function of all of these. We split 
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Figure 10. Stacked SPIRE luminosities as a function of g — r colour, rcdshift and absolute magnitude M r . Layout as in Fig. [8] Error bars 
are the statistical la errors in the bins as described in Section l3.3l Fluxes in the red bin have been corrected for the lensing contribution 
as described in Section 14.21 and error bars include the associated uncertainty. The thick grey line is the same from left to right, and is 
the linear least-squares fit to the results for the lowest-redshift blue galaxies. 



the sample by colour first, then by redshift and finally divide 
into mass or magnitude bins, but each bin can still contain 
a relatively broad range of redshifts. Within each bin there 
will be a strong degeneracy between redshift and M r , sim- 
ply because M r is a strong function of redshift. In Fig. [T3] 
we plot stellar mass against M r with the points colour-coded 
by redshift, showing that redshifts increase steadily from left 
to right, with decreasing M r . A narrow range in M r would 
select a narrow range of redshifts, while a similarly narrow 
range in mass selects a much broader range of redshifts. 



The effect of this on the blue bins in Figures [10] and QT] 
is a tendency for the data points of different redshift bins to 
lie along the same relation of L250 as a function of M r . The 
degeneracy is (partially) broken when splitting by stellar 
mass, thus separating out the trends with redshift and with 
mass in Fig. 1121 This effect is much less noticeable in the red 
bin simply because the redshift evolution is much stronger 
while the mass dependency is very weak in the red sample. 
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Figure 11. Stacked 250 fim luminosity as a function of NUV — r colour, redshift and M r . Error bars are the statistical la errors in the 
bins as described in Section 13.31 Data and errors in the red bin incorporate the correction for lensing. The thick grey line is the same 
from left to right, and is the linear least-squares fit to the results for the lowest-redshift blue galaxies. 
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Figure 13. Stellar masses of the sample as a function of abso- 
lute magnitude and coloured by redshift, showing the degeneracy 
between M r and z resulting from an r-band selection. To guide 
the eye, the slope of the grey line indicates direct proportional- 
ity between mass and luminosity, i.e. log 10 M s t a r = — 0.4M r + C 
(for this line C = 1.0). The spread in the data perpendicular to 
this line reveals the broad range of mass-to-light ratios which is 
responsible for the differences between stacking by M r and by 



5 DISCUSSION 

5.1 Dust temperatures and SED fitting 

In Section 14.41 we stated that assumptions about the dust 
temperature and had negligible effect on the fc-corrections 
to SPIRE fluxes at these low redshifts. Hence we chose 
to use the same SED model to compute monochromatic 
luminosities, assuming a single-component greybody with 



Tkust = 18.5 K and /3 = 2.0. However if we want to infer 
properties of the full IR SED these considerations are much 
more important. 

We fitted single-component SEDs with /3 = 2.0 to the 
stacked SPIRE fluxes in each bin, shifting the observed 
wavelengths by (1 + z) using the median redshift in the cor- 
responding bin. Fluxes in the red bin were first corrected 
for the predicted lensing contamination as described in Sec- 
tion 14.21 The effect of this is to increase the fitted tempera- 
tures in the red bin by around 1 — 3K, which is small com- 
pared with the range of temperatures observed, although 
the errors on temperatures are significantly increased. The 
fitting was carried out using the IDL routine mpfitfunU^I 
which performs Levenberg-Marquardt least-squares fitting 
to a general function. Best-fit values of the free parameters 
(temperature, normalisation) are returned with formal la- 
errors computed from the covariance matrix. Some exam- 
ples of the fits are plotted in Appendix [Cj showing a range 
of fitted temperatures. The derived temperatures depend on 
the assumption of a fixed emissivity parameter (/}). Varying 
this as a function of optical colour and/or stellar mass could 
to some extent account for the variation in submm colours, 
which we interpret as a temperature variation. However the 
variation in /3 would need to be severe (A/3 > 1) to fully 
account for the trends in the stacked colours in Fig. [9] It 
therefore seems likely that the cause for these variations is 
either temperature alone or a combination of /3 and temper- 
ature. 

Figure [14] shows the results of all the temperature fits 
as a function of colour, stellar mass and redshift. There are 
strong deviations in some bins from the value that we have 
been using, and these show strong dependence on colour and 
stellar mass. In the blue bin we see that dust temperature 
is tightly correlated with stellar mass in all redshift bins, 



10 mpfitfun available from Craig Markwardt's IDL library: 
http: / /cow. physics . wise . edu/~ craigm/idl/idl .html 
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Figure 12. Stacked SPIRE luminosities as a function of g — r colour, redshift and stellar mass. Error bars are the statistical la errors 
in the bins as described in Section 13,31 Data and errors in the red bin incorporate the correction for lensing. The thick grey line is the 
same from left to right, and is the linear least-squares fit to the results for the lowcst-redshift blue galaxies. 



with a peak at around 6 x 10 Mq, but galaxies of higher 
masses appear to have colder dust. The temperature distri- 
bution in the the green bin is less well correlated and very 
noisy, but again there is evidence that the warmest galax- 
ies are towards the middle of the mass range. There is no 
clear evolution with redshift. As with the luminosity results, 
the red bin appears very different from the other two, and 
there is no evidence for the temperature to increase with 
stellar mass. The most important result would seem to be 
that temperatures are generally much lower in the red than 
the blue bins. Over the mass range in which the bins overlap 



(3 X 10 9 < M st ar < 3 x 10 11 ), the mean (standard deviation) 
of the temperatures in the blue bins is 22.7 (2.9) K, com- 
pared with 19.4 (2.8) K in the green and 16.1 (2.9) K in the 
red. The difference in the means is statistically significant, 
since an unpaired i-test gives a probability of 10 -12 that 
the means of the red and blue bins are the same. The t-test 
assumes that the errors on each of the measurements are 
independent, but this is not true since the errors on the red 
stacks are dominated by the error on the lensing correction. 
The mean temperature error of the red stacks is 3.1 K, which 
means that including the lensing uncertainty, the blue and 
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red mean temperatures are only different at the 2.1a level. 
The incidence of colder dust in redder galaxies may be ex- 
plained by the relationship between dust temperature and 
the intensity of the interstellar radiation field (ISRF). The 
colour temperature of the ISRF is directly related to the 
stellar population. Old stars produce less UV flux, and so 
heat the dust less, which causes galaxies dominated by older 
stars to have cooler dust. Such a temperature differential is 
therefore consistent with the notion of red galaxies being 
passive. Meanwhile, colder dust temperatures in the least 
massive blue galaxies is consistent with these galaxies having 
more extended dust disks in comparison with their stellar 
disks, since the ISRF becomes weaker at greater galacto- 
centric radii. Evidence for an extended dust disk in at least 
one lo w mass system has been reported bv lHolwerda et ail 
although larger samples would be needed to judge 
whether this is a widespread phenomenon. 

It is possible when fitting the SPIRE bands that the 
SED shape could be biased by differential effects between 
the three bands. In particular the 500 /im band is the 
most affected by confusion and blending, as well as being 
the noisiest, and is also potentially subject to contamina- 
tion from other emission mechanisms, including synchrotron 
from within the galaxies, and extended radiation from the 
Sunyaev-Zel'dovich effect in galaxy clusters (although we 
note the latter should have been removed in background 
subtraction). To check for such bias, we also tried fitting 
only the 250 and 350 /im data, and found that the derived 
temperatures and trends were not significantly different. 



5.1.1 Bolometric luminosities 

It can be useful to consider the total IR luminosities Ltir, 
(8 — 1000 /im) of galaxies as this allows some comparison 
between observations at different IR wavelengths and be- 
tween data and the predictions of models. In all cases this 
entails making assumptions about the shape of the SED, 
which must be interpolated - and indeed extrapolated - 
from the limited photometric data available. In this par- 
ticular case we are limited to just three SPIRE bands, all of 
which lie longward of the peak in the SED and as such do 
little to constrain the warmer end of the SED at A < 100 /im. 
This is why they are well fitted by single-component SEDs, 
representing a single component of cold dust. In contrast, 
the TIR luminosity is highly sensitive to emission from the 
hotter components of dust, especially the > 30 K dust as- 
sociated with Hn regions, which is heated by UV radiation 
from hot young stars. 

Bearing in mind these limitations, we nevertheless con- 
sider it useful to make some attempt at estimating the TIR 
luminosities representative of our stacked samples. Since our 
sample is thought to be dominated by normal star-forming 
and quiescent galaxies, we need to choose an appropriate IR 
SED template. A com monly us ed set of templates is that of 
I Chary fc Elba3 (|200ll . hereafter ICEOll '). These templates are 
based on libraries of mid- and far-IR templates represent ing 
a range of SED types (from normal spirals to ULIRGa 11 !) 
fitted to data on ~ 100 local galaxies at 6.7, 12, 15, 25, 



11 Ultraluminous IR galaxies; Ltir, > 10 12 Lq. 



60, 100 and 850 /mi[3 From these we select the most ap- 
propriate template for each stack by computing chi-squared 
between each of the templates and our rest-frame (lensing 
corrected) SPIRE luminosities, and assign to each stack the 
Ltir of the template with the minimum chi-squared. Results 
are shown in Fig. ll5f a). Errors on Ltir were estimated with 
Monte-Carlo simulations using the la errors on the SPIRE 
luminosities and re-fitting the templates 200 times to obtain 
the la e rror ba r on the template Ltir. 

The ICEOll templates are of limited value for our sam- 
ple because they are fitted to IRAS and SCUBA data for a 
relatively small sample of local galaxies. The necessity for 
IRAS detections means that the galaxies in their sample 
may have been biased towards hotter SEDs, and may not 
be representative of the larger population sampled in this 
work. As an alternative we can compare the results of us- 
ing the lCEOll templates with a set of te mplates modelled o n 
the H- ATLAS SDP source catalogue (|Smith et all 1201 lbl ). 
There is a danger that the opposite bias is active here, 
since the templates are based on sources selected at 250 /im, 
which are more likely to have cold S EDs. However by com- 
paring the H- ATLAS L250 LF from iDllI with the range of 
L250 of optical galaxies (Figures I10H12[) we see that the lu- 
minosity ranges spanned by the two surveys are remark- 
ably similar, implying that the SEDs of H- ATLAS sources 
could provide a reasonable representation of an optically se- 
lected sample. We use a single tem plate based on the m ean 
of all H- ATLAS SED models from lSmith et all (|2011bl L In 
Fig. I15f b) we show the results of fitting this template to 
our stacked SPIRE luminosities, minimising chi-squared to 
obtain the correct normalisation and integrating the SED 
from 8-1000 /an to obtain Ltir (errors were estimated from 
Monte-Carlo simulations using the la errors o n the SPIRE 
luminosities in the same way as for the ICEOll templates) . 

The results of the two se ts of templates are strik- 
ingly different, with the ICEOll models suggesting signifi- 
cantly higher luminosities, reaching the level of 'luminous IR 
galaxies' (LIRGs; Ltir > lO n L ) at z > 0.22 or M star > 
2 x 10 10 M Q . The H- ATLAS templates are much colder so 
give much more moderate luminosities, with around five 
times lower bolometric luminosity for the same L2so- With 
only the SPIRE data to constrain the SED we cannot con- 
clusively say that either set of templates is better suited 
to describing the optical sample, although for the reasons 
outlined above we believe that the H-ATLAS templates are 
more likely to be suitable. The addition of data points at 
shorter wavelengths, from PACS in the FIR and WISE in 
the MIR, would permit a much more accurate derivation of 
the bolometric luminosity; we leave this for a future study. 

Both parts of Fig. [TS] show trends in the bolometric 
luminosities that are similar to those seen in the monochro- 
matic SPIRE luminosities. That is to say there is a clear 
evolution with redshift and that this is much stronger in the 
red than the blue sample. We can quantify this evolution in 
the form L(z) oc (l+z) a for the mass range in which the bins 
overlap. We fit this function in log-space by chi-squared min- 
imisation, using the IDL routine linfit to fit the bins which 
fall in the mass range AW = 2 - 7 x 10 10 M . We find that 

12 ICEOll templates were obtained from 

|http : //www . its . caltech. edu/- rchary/ 



© 0000 RAS, MNRAS 000, 000-000 



22 N. Bourne et al. 



BLUE: O.OOS g-r<0.05-0.026M r 



GREEN: 0.05-0.026M.S q-r<0. 1 5-0.026M, 



RED: 0.15-0.026M r S g-r<1.00 



30.0 



25.0 



20.0 



15.0 



10.0 




X 0.01 s 


z 


<0.12 


O 0.12S 


z 


<0.17 


A0.17S 


z 


<0.22 


□ 0.22S 


z 


<0.28 


0.28S 


z 


<0.35 



9 10 11 

l°9,o(M stor [MJ) 



9 10 1 

l°9,o(M stor [MJ) 



9 10 1 

l°g 10 (M slor [MJ) 



Figure 14. Results of fitting single-component greybodies with /? = 2.0 to the observed (not k— corrected) fluxes in each bin to estimate 
dust temperatures. Fluxes in the red bin were corrected for lensing before fitting. Error bars are the la errors on fitted temperatures 
computed by mpfitfun. 



Ltir from the H-ATLAS SEDs evolves with a = 4.1±0.2 in 
the blue sample; a = 1.3±0.3 in the green; and a = 6.9±1.0 
in the red. The same evolution is also seen in L2so{z) in 
the same mass bins: this is described by a — 4.0 ± 0.2 
(blue); a = 1.1 ± 0.4 (green); and a = 7.7 ± 1.6 (red). It 
appears counter-intuitive that the intermediate green bin 
should show the least evolution. The likely reason for this is 
the aforementioned possibility for the green sample to probe 
different populations at different redshifts (see Sections 14.51 
and !4.1|l . Any sign of genuine luminosity evolution would be 
counteracted by sampling a less luminous population (e.g. 
with more passive red sequence galaxies in the green bin) at 
higher redshifts. We must however note that any evolution 
in the submm SED of any of these samples (blue, green or 
red) would render these single template fits unreliable. 

The evo lution in Lt i r of n ormal galaxies was also 
observed by lOliver et al. | (|2010bh for a large sample se- 
lected in the optical and NIR with redshifts between 
— 2. Dividing their sample by redshift, stellar mass and 
op tical class (each derived fro m broadband SED fitting 
of H owan-Robinson et al.l 120081 ) they stacked into 70 and 
160 (im Spitzer images and showed that both 'blue' galaxies 
(with spiral-like SEDs) and 'red' galaxies (with elliptical-like 
SEDs) increased in specific IR luminosity (i.e. LTiR/M star ) 
as a function of redshift. The evolution in specific IR lu- 
minosity of all galaxies in their sample increased as (1 + 
^4.4±o,3 (independent of stellar mass), which is nearly iden- 
tical to our result for blue galaxies. When they split the 
sample into red and blue colours, they also found that 
red galaxies evolved more strongly, but only with the in- 
dex 5.7 ± 2.5 compared with our 6.9 ± 1.0. Their blue sub- 
sample evolved with the index 3.4 ± 0.3 compared with our 
4.1 ± 0.2. The agreement is not exact but general trends 
with redshift and colour are certainly compatible between 
the two samples. Assu ming a correlatio n between IR lumi- 
nosity and SFR (e.g. iKennicuttl 1 19981 ). we can also draw 
paral l els with other ffjpfee r-stacking (such as lMagnelli et all 
120091 ; iDamen et alj|2009al lbh as well as radio-stacking stud- 



ies dDunne et all l2009al ; IPannella etafl 120091 ; iKarim et all 
l201l| ). all of which have shown similar dependence of (spe- 
cific) SFR on stellar mass and redshift in NIR-selected mas- 
sive galaxies covering larger redshift ranges (up to z ~ 3). 
These studies have variously reported redshift evolution in 
specific SFR with indices (a) ranging from 3.4 to 5.0, all 
comparable with the luminosity evolution of our blue sam- 
ple. It is unsurprising that our blue sample generally agrees 
with other samples selected by rest-frame optical light with 
no regard to colour, since our blue selection is by far the 
largest of our colour bins, comprising 50 per cent of our 
sample. 

It is well reported in the literature that there is 
strong evolution in the IR LF out to at least z = 1 



dSaunders et all 19901; 


Blain et all 


199S 




Pozzi et all 2004; 


Le Floc'h et al. 2005; 


Eales et all 


2009 


2010al; Dve et all 


20ld; GruDDioni et all 20101; Rodiehiero et all |2010|; iDllI; 


Goto et al.l201ll; Sedgwick et al. 2011). 


This requires an in- 



crease in the luminosity of the brightest (> L\ K ) galaxies, 
leading to an increase in the numbers of LIRGs and ULIRGs 
at higher redshifts. Our results show that this evolution at 
low redshifts also occurs in ordinary galaxies well below the 
LIRG threshold; these are the galaxies that dominate the 
number density. Such an evolution in the IR luminosities of 
all galaxies leads naturally to an evolution in the character- 
istic luminosity L* . Exactly this sort of evolution in normal 
(i.e. non-merging) star- forming galaxi e s is p redicted by the 
semi-analytic model of iHopkins et al] l|2010t ). essentially as 
a result of an e volving gas fracti on and using the Schmidt- 
Kennicutt law (|Kennicuttl Il99cf l. IDllI also show that an 
evolving gas fraction is required to explain the luminosity 
evolution in the H-ATLAS sample, based on the chemical 
evolution model of Gomez, H.L. et al. (in preparation). 



5.2 The cosmic spectral energy distribution 

Having discussed evolution in the IR luminosity density of 
the Universe, it is natural to consider the local luminosity 



© 0000 RAS, MNRAS 000, 000-000 



H-ATLAS/GAMA: dust in optically selected galaxies 23 



BLUE: O.OOS g-r<0.05-0.026M r 



(a) CEOl templates 
GREEN: 0.05-0.026M r S g-r<0. 1 5-0.026M, 



RED: 0.15-0.026M r S g-r<1.00 




9 10 11 

log,o(M stor [MJ) 



9 10 

l°9,o(M stor [Mj) 



9 10 1 

l°9,o(M slor [MJ) 



1 1 



10 



BLUE: O.OOS g-r<0.05-0.026M r 



X 0.01 s 


z 


<0. 1 2 


O 0.12S 


z 


<0.17 


AO. 175 


z 


<0.22 


□ 0.22S 


z 


<0.28 


0.28S 


z 


<0.35 



4* 



9 10 11 

l°9io(M star [MJ) 



(b) H-ATLAS templates 
GREEN: 0.05-0.026M r S g-r<0. 1 5-0.026M, 



9 10 11 

log,o(M star [MJ) 



RED: 0.15-0.026M r S g-r<1.00 



A- 



9 10 11 

log,o(M star [MJ) 



Figure 15. Integrated L TIR (8-1000 /im) of (a) ICharv fc Elbaj | |2001|) and (b) H-ATLAS jSmith et alj|2011bh templates fitted to the 
stacked SPIRE luminosities in each bin as described in the text. Error bars are la errors estimated from Monte-Carlo simulations using 
the lcr errors on the SPIRE luminosities. Fluxes in the red bin have been corrected for the lensing contribution as described in Section l4.2l 
and error bars include the associated uncertainty. Note that in panel (a) the first point in the blue sample (i.e. lowest mass, lowest redshift 
bin) had SPIRE luminosities that fall below all of the CEOl templates, and so the luminosity of the faintest template is used as an upper 
limit. 



density (at 2 = 0), since this provides a reference point for 
similar measurements at higher redshifts. In Section [4.31 we 
calculated the integrated intensity of low redshift galaxies 
and showed that they contribute a small fraction of the CIB 
at submm wavelengths. Building on this result, we can esti- 
mate the z = cosmic SED at submm wavelengths, i.e. the 
integrated luminosity of all galaxies at z = 0. To do this we 
make use of the completeness-corrected integrated intensi- 
ties in the range 0.01 < z < 0.12 in Tabled but apply k- and 
e-corrections to account for the redshifted wavelengths and 
luminosity evolution respectively. We divide by the comov- 
ing volume of the redshift bin (V c ) to obtain the luminosity 



per unit comoving volume: 

jO K(z) 4-7T . . 

1 + 2 Vc 

where I v is in Wm -2 sr _1 Hz -1 , d_t is in m, V c is in 
Mpc 3 , and the luminosity vL v is expressed in units of 
WMpc" 3 /i 70 (hro = H o /70 km s" 1 Mpc" 1 ). Note that while 
V c = 0.523 Gpc 3 represents the total comoving volume of the 
bin, d,L — 382.5 Mpc is the luminosity distance of the me- 
dian redshift in the bin, (z) = 0.084. We use the median 
values oi K{z)/(l + z) for this redshift bin: 0.8094 (250 /im), 
0.7586 (350 /im), and 0.7259 (500 /im); and for the evolu- 
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tion at all submm wavelengths we take the fitted function 
£250 {z) oc (l + z) a from Section r5.ll assuming that the shape 
of the rest-frame SED does not evolve with redshift (which 
is supported by our non-evolving temperature results). We 
assume a = 4.0 ±0.2, the value derived for the blue sample, 
since these represent roughly half of all galaxies; the evolu- 
tion for all galaxies may be slightly stronger (red galaxies 
seem to evolve more strongly, although the green sample 
evolve less) but an index around 4 is consistent with other 
results in the literature (see discussion in Section f5 . 1 [) . The 
evolution from z — — 0.084 (the median redshift of the 
bin) is therefore a factor (1 + 0.084) 40 = 1.38, hence the 
correction is e(z) = 0.72. 

We thus calculate the luminosities of the cosmic SED 
at z = to be 3.9 ± 0.3 x 10 33 , 1.5 ± 0.1 x 10 33 and 
4.3 ± 0.3 x 10 32 WMpc~ 3 /t 7 o at 250, 350 and 500 /mr respec- 
tively. Errors are dominated by the 7% calibration error on 
the SPIRE fluxes, which is correlated across the three bands. 
These results agree very closely with the submm luminosities 
predicted from GAMA data by Driver, S. et al. (in prepara- 
tion), by calculating the total energy absorbed by dust in the 
UV-NIR and ass uming it is reprocesse d as FIR emission with 
templates from I Dale fc Heloul (|2002h . They are also close 
to the pre-Herschel-era prediction of ISerieant fc Harrison! 
|2005l ). based on modelling the SEDs of IRAS sources with 
SCUBA submm measurements. Using equation (7) of that 
paper, the predicted luminosities at 250, 350 and 500 /im 
are 4.52 x 10 33 , 1.43 x 10 33 and 3.48 x 10 32 W Mpc" 3 h 70 
respectively. Our measurements are within 3a of these val- 
ues, although they arguably suggest that the slope of the 
cosmic SED may be a little shallower than the prediction. 
These measurements are independent of the SEDs and tem- 
peratures assumed (since fc-corrections are small) and of the 
lcnsing assumptions (the lensed flux is negligible at z < 0.1). 



5.3 Evolution of dust masses 

Dust mass is a quantity which we can expect to constrain 
much more accurately than Ltir using the SPIRE luminosi- 
ties, because the cold component that they trace is tho ught 
to dominate the total dust mass (IDunne fc Ealesll200ll , and 
references therein). We therefore estimate the mass of the 
cold dust component described by our greybody fits and use 
this as a proxy for the total dust mass, assuming that any 
warmer components have a negligible contribution to the 
mass. 

The dust mass as a function of temperature Td us t is 
estimated from the 250 /im flux using equation ([9|: 



Mdust = 



5250 Dl K(z) 



K250 B(V250, Tdust) (1 + z) 



(9) 



We use a dust mas s ab sorption coefficient at 250 fim of 
K250 = 0.89 m 2 kg -1 (|Plll . and references therein). 

Dust mass results depend strongly on the temperatures 
assumed (although not as strongly as the bolometric lumi- 
nosities). They are therefore subject to our assumption of 
constant emissivity index (/?), as well as the assumption of 
a constant absorption coefficient (k). Any variation of k with 
redshift, stellar mass, optical colour or indeed with dust tem- 
perature or dust mass itself could alter the trends that we 
see in dust mass. 

In Fig. [16] we show the dust masses derived using the 



fitted temperatures from Fig. 1141 The dust mass is seen 
to range from around 2 x 10 6 to 8 x 1O 7 M0 across the 
sample. In all bins dust and stellar mass are correlated, 
but this weakens slightly with increasing redshift and/or 
stellar mass. There is a definite evolution towards higher 
dust masses with increasing redshift. Following the evolu- 
tionary form Mdust(z) oc (1 + z) a we fit data in the range 
Mstai = 2 - 7 x 10 10 M with slopes of a = 3.9 ± 1.7 for the 
blue sample; a = 3.0 ± 3.3 for the green; and a = 6.8 ± 4.6 
for the red. The slopes of the evolution are thus similar to 
the evolution in luminosities (as might be expected from the 
lack of evolution in temperature). Given the error bars we 
cannot claim to detect evolution in the dust masses of red 
galaxies, despite having a significant detection of evolution 
of their luminosities. The reason for this is that the uncer- 
tainty of the lensing contribution increases the uncertainty 
in the fitted temperatures. If we assume that dust temper- 
atures of the red sample do not evolve (as they do not for 
the blue and green samples) then the evolving luminosities 
must result from dust mass evolution. However this may not 
be a valid assumption if the composition of the red sample 
changes with redshift. 

If we ignore the highest redshift red bin, which has the 
largest errors due to lensing, then the results seem to sug- 
gest that the difference in dust mass between the red and 
blue galaxies is weaker than the difference in luminosity, for 
a given stellar mass and redshift. The dependence of lumi- 
nosity on colour therefore appears to be driven more by the 
temperature than the mass of the dust. We note that if the 
same temperatures were used in deriving the dust mass in 
every bin then the dust masses would be directly propor- 
tional to 1/250 and would follow the trends seen in Fig. 1121 
This shows the vital importance of having photometry at 
multiple points along the SED, without which it would be 
impossible to constrain the SED shape and inferences about 
dust mass evolution would have to assume a constant tem- 
perature. 

Significantly, this analysis suggests that the cold dust 
masses of red and blue galaxies are not strikingly discrepant 
in the stellar mass range ~lx 10 1() - 2 x IO^Mq. Previ- 
ous studies fitting two-component dust models to normal 
spiral galaxies in the local Universe have derived ranges of 
cold dust masses com parable to our sa mple: 3 x 10 5 - 1 x 10 s 
jPopescu et alj|2002h ; 2 x 10 7 - 1 x 10 9 dStevens et alj|2005h ; 
and 4 x 10 6 - 6 x 10 7 M Q (|Vlahakis et ajj|2005h . Dust masses 
measured from fits to the cold dust in local ellipticals (which 
should be akin to o ur red sample) are often a little lower: 
2 x 10 5 - 2 x 10 6 (|Leeuw et al.l 120041); 4 x 10 4 - 5 x 10 7 



dTemi et alj|2004h: 2xl0 4 -2xl0 7 (ISavov et al.ll2009h . How- 



ever, both IVlahakis et al.l (|2005l ) and IStickel et all (|2007h 
reported little or no significant difference in the typical 
dust masses of galaxies of different Hubble types (includ- 
ing spheroidals, spirals and irregulars), although Stickel et 
al. reported that spheroidal and irregular types reached sig- 
nificantly lower dust masses. One particular issue noted by 
Vlahakis et al. was the possibility of contamination of their 
850 /im fluxes with synchrotron emission, which would lead 
them to overestimate a few of their dust masses. We note 
the caveat that in contrast to our unbiased sample, the ref- 
erences in this paragraph were studies of individual galaxies 
selected variously with IRAS or ISOPHOT in the FIR, or 
SCUBA in the submm, so the range of dust masses sampled 
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would not have been comp lete (with the excep tion of the 
optically selected sample of lPopescu et al]|2002T l. 

iPopescu et all (|2002r i observed colder dust in later Hub- 
ble types (which might be expected to be the bluest galax- 
ies). Their sample of spirals would probably reside entirely 
within our blue bin so such a trend would not be apparent 
between our colour bins if it does not extend beyond the 
blue cloud. However there is a correlation between Hubble 
type and mass; later types have lower stellar masses, so our 
observation of a strong correlation between stellar mass and 
dust temperature in blue galaxies is entirely consistent with 
the findings of Popescu et al. 

M eanwhile, results fro m the Herschel Reference Survey 
(HRS: iBoselli et alllioiOrl ) indicate that early type galax- 
ies (E+S0+S0a) detected by Herschel in a volume-limited 
sample of the local Un iverse have dust ma sses in the range 

1 x 10 5 - 2 x 10 7 M (jSmith et al]|2011ch - although they 
only detected 34 per cent of ellipticals and 61 per cent of 
SO's. Their sample have NUV — r colours that place them 
in our 'red' bin, yet their dust masses are much lower than 
the typical dust masses that we find for red galaxies. This 
is perhaps due in part to the higher derived temperatures 
of the HRS sample: Td ust = 16 — 32 K; with a mean of 24 K 
in comparison with our mean of 16.1 K for the red sample 
(both assuming /3 — 2). Smith et al. concluded that the dust 
masses of SO's were around 10 times lower than those of the 
HRS spirals, while those of ellipticals were 10 times lower 
again (for the same stellar mass), which seems to contrast 
with our results. 

iRowlands et al.1 (|2012T ) studied early type galaxies de- 
tected in H- ATLAS and found dust masses mostly between 

2 x 10 7 and 2 x 10 s M , with a mean dust mass similar 
to that of spirals (5.5 x 10 7 Mq). However, the stellar mass 
distributions of spirals and e arly types were v ery different. 
The NUV — r colours of the IRowlands et al.1 spiral sample 
lie mostly within our blue bin, and the early types mostly 
within our green bin. Their redshift range is similar to ours 
(z < 0.3), and the mean dust and stellar masses of their 
spirals/early types lie within the locus of our blue/green 
bins in Fig. 1161 However, derived dust masses depend on the 
temperature ass umed. The cold dust temperatures fitted by 
IRowlands et al.l range from 15 to 25 K, while the temper- 
atures in our blue/green samples for the same stellar mass 
range (> 1O 1O M ) are between 15 and 28 K. The correspon- 
dence is not exact but the ranges are similar so average dust 
mass results should be comparable. For reference, chang- 
ing the temperature from 15 to 25 K results in a drop in 
the derived dust mass by a factor 5, which is similar to the 
range of dust masses across the redshift range in Fig. 1161 
IRowlands et al.1 (|2012T l compared their H-ATLAS-detected 
early types with a control sample of undetected early-types 
selected to have a matching distribution of redshifts and r- 
band magnitudes. They concluded that the detected early 
types were unusually dusty compared with the control sam- 
ple, and could be undergoing a transition from the blue cloud 
to the red sequence. It seems likely that those objects do in- 
deed comprise the top end of the dust mass distribution 
(at a given stellar mass and redshift), although they do not 
appear to be exceptional outliers when compared with the 
median dust mas ses in our sampl e. This is surprising when 
we consider that IRowlands et al.l found the detected early 
types to have around 10 times as much dust as typical early 



types: why are they not also outliers compared to typical red 
galaxies? The answer may be that typical red galaxies are 
generally dustier than typical early types, which supports 
the notion (as discussed in earlier sections) that our red sam- 
ple is comprised of a mixture of different populations. The 
red bin is likely to contain most of the early type galaxies in 
the sample volume, and if these are relatively dust-poor then 
there must be a substantial population of red dusty galaxies 
boosting the median dust masses in the red bin. This could 
also explain the d iscrepancy between our red sample and the 
HRS early types jSmith et al.ll2011d\ 

There is also the possibility of an environmental factor 
in the offset between the HRS results and our own. Many of 
the galaxies in the HRS sample reside in the Virgo cluster, 
while most of the galaxies in GAMA will be in lower density 
environments. The lower dust masses of HRS early types 
compared with GAMA red galaxies could be due to early 
types in clusters being dominated by passive red-sequence 
systems, while red galaxies in lower density environments are 
more likely to be dusty. Such a division is indeed suggested 
by the higher detection rate of early types outside of V irgo 
in HRS, compared to those inside (|Smith et al.ll201ld ). It 
is however unclear whether such an effect could be strong 
enough to fully explain the discrepancy we find. 

On the other hand, if the lensing contamination is 
slightly greater than we have predicted, then the derived 
dust masses of our red sample could be a lot lower; how- 
ever this is only likely to affect the higher redshift bins 
due to the wea k lensing efficiency at lower redshifts. The 
IRowlands et al.l sample is less likely to be biased by strong 
lensing than our sample is, because they excluded detec- 
tions. The HRS results are unlikely to be biased by lensing 
because their fluxes were much higher than the red galaxies 
in our sample, and the low dust masses and high tempera- 
tures they derive (relative to our own) argue against their 
results being biased by lensing. 

5.3.1 Dust-to-stellar mass ratios 

In Fig. [17] we plot the ratio of dust to stellar mass across 
the sample, which shows several interesting features. Firstly 
there is in general a strong correlation with stellar mass: 
the more massive galaxies have smaller dust-to-stellar mass 
ratios. The correlation appears steeper for red galaxies of 
the highest masses in each redshift bin, but this is not so 
obvious in the blue and green samples which do not reach 
to quite such high stellar masses. Nevertheless it seems not 
unreasonable to observe that the most massive red galaxies, 
many of which will be passively evolving giant ellipticals, 
have especially low dust-to-stellar mass ratios. 

It is worth pointing out the potential for a negative cor- 
relation between M duBt / M Bt!l T and M ata r to be produced arti- 
ficially in binned data. This can occur if there is a large range 
of stellar masses with large errors, even when there is no cor- 
relation between Md UB t and M a t ar , since a bin that selects 
data with low M star also selects those with high Md us t /A/star • 
In general the slope of the measured correlation could be af- 
fected by this artificial phenomenon; however we can be sure 
in this case that the trends are real because they can also 
be discerned in the median Md us t values in Fig. 1161 and in 
any case our stella r mass errors are small (AlogM sta r ~ 0.1; 
iTavlor et ai]|201ll l. 



© 0000 RAS, MNRAS 000, 000-000 



26 N. Bourne et al. 



The aforementioned redshift evolution is very apparent 
in Fig. 1171 and although we naturally select higher stellar 
masses at higher redshift, the dust masses in the sample rise 
more rapidly resulting in an increasing dust-to-stellar mass 
ratio with redshift. Using the (1 + z) a model once again 
we find that the evolution is consistent with the slopes de- 
rived for the dust mass evolution. This evolution in dust 
mass echoes the results of Dll, who found a strongly evolv- 
ing dust mass function (DMF) in H- ATLAS sources up to 
z ~ 0.4, a s well as results from oth er surveys reaching higher 
redshifts (lEales et al.ll2009l,l2010al) Using dust masses from 
SED fitting bv ISmith et all (|2011bl ). [Dill showed that the 
characteristic dust mass (M| ust ) of the H-ATLAS DMF in- 
creases from 3.8 x 10 7 M Q at z < 0.1 to around 2.1 x 10 s M 
at z ~ 0.35 (although they note that this does not measure 
the true evolution because there is an accompanying fall in 
the characteristic density cjf). This range is indicated by 
the shaded region in Fig. 1161 In all of the bins, our typical 
dust masses reach lower than the mini mum mass sampled 
in equivalent redshift slices in the |P 111 DMF, but the fact 
that we see evolution indicates that galaxies of a given stel- 
lar mass shift up the DMF at increasing redshifts. We see 
this happening in galaxies of all colours and stellar masses, 
indicating that the evolution in the DMF is the result of 
changing dust masses in all galaxies, both passive and star- 
forming. The evolutio n (ar ound a factor 3 — 4) we see is 
similar to that seen bv lDlll . who fitted two-component dust 
masses temperatures using a detailed physically-motivated 
SED model. 

It is also interesting to compare Md ua t/Af a tar in our 
results with the detected H-ATLAS galaxies in iDlll . The 
H-ATLAS galaxies were found to have higher dust-to-stellar 
mass ratios than predicted by models, ranging from 2 x 10~ 3 
at z < 0.1 to 7 x 10~ 3 at z ~ 0.3 (this range is shaded in 
Fig. I17|l . Our results are typically lower but the strong de- 
pendence on stellar mass means that they span a very wide 
range; many of the blue galaxies in the lower mass bins have 
much higher dust/stellar mass ratios than the H-ATLAS 
sample for the same redshifts. It is perhaps not surpris- 
ing that we sample a much wider range than the H-ATLAS 
sources since our sample selection criteria are independent 
of dust content. The dust /stellar mass results have implica- 
tions for understanding the dust production mechanism, as 
we can show by comparing the results from a chemical evolu- 
tion model with the parameters obtained for the H-ATLAS 
galaxies (Gomez, H. L. et al. in preparation). The m odels 
(based on the framework in lMorgan fc Edmu nds 2003) show 
that values of Md US t/Af s tar > 10 -3 cannot be achieved with 
a purely stellar source of dust . Even including dust produc- 
tion in supernovae ejecta (e.g. iRho et al . 2008: Dunne et al.1 



l2009bl : iGomez et al.ll201ll : iMatsuura et al.ll201ll ). models re- 
quire the condensation efficiency in the ejecta to be close to 
100 per cent to reach the high values of Md us t /A/star ~ 10~ 2 
seen both in the H-ATLAS detected sample and in the 
lowest-mas s blu e galaxies in our sample. However, as dis- 
cussed by ID 111 and Gomez et al. (in preparation), such 
high dust yields from supernovae are difficult to produce, 
leading to the invocation of alter native explanations such 
as dust grain growth in the ISM jDraine fc Sarpeterlll979l ; 
ID wek fc Scaioll 19801; iDrainel 1 1 990l; fedmund j|200 if) or a to 
heavy IMF (e.g. lHaravama. Eisenhauer fc Martina l200i 
The models also indicate that the low mass galaxies with 



high Mdu S t/M s tar are less efficient at turning their gas into 
stars (compared to high mass sources with low Md us t/Af s tar). 
In other words, low mass systems have a longer star forma- 
tion timescale, so although they produce less dust mass per 
year from stars, there are more metals and dust in the ISM 
for longer, in comparison to massive galaxies (Gomez et al. 
in preparation). 

Moving on to higher redshifts, ISantini et al.1 (|2010h 
showed that the dust/stellar mass ratios of ordinary low 
redshift galaxies are much lower than those of high-redshift 
submm galaxies (SMGs) from the Herschel-PEP survey. Our 
Afdust /M B tar values from stacking are consistent with their 
sample of low-redshift spirals from the SINGS survey (with 
typical stellar masses of - 10 n M fl ) Their dus t masses were 
derived by fitting GRASIL models (|Silva et al.ll 19981 ') to pho- 
tometry spanning the FIR/submm SED, and are consistent 
with the value of /3 = 2 that we ha ve assumed. O ur results 
therefore support the conclusion of ISantini et al.l that high- 
redshift SMGs have much higher dust content (by a factor 
~ 30) than local spiral galaxies. This offset is also consistent 
wi th the comparison between low-redshift H-A TLAS sources 
in ID 111 and high-redshift SCUBA SMGs in iDunne et af] 
|2003l ). It has now become clear that the SMGs detected 
in early submm surveys are exceptionally dusty systems in 
comparison with the low redshift galaxy population. 

There is evidence th at the dust /stellar mass ratio corre- 
lates w ith specific SFR (|da Cunha et alJl201Ct ISmith et al.l 
l2011bl . Gomez, H. L. et al. in preparation), so these results 
imply that the least massive galaxies at low redshift are 
the most actively star-forming, and that all galaxies become 
more active towards higher redshifts. At z ~ 0.3, galaxies 
with 2 x 10 10 Mq of stars have the same JV/d us t /M s tar as 
galaxies a quarter of that size do at z < 0.1. This is con- 
sistent with the picture of downsizing, in which the specific 
star formation rates of high mass galaxies pe ak earlier in 
the U niverse than those of low-mass galaxies (|Cowie et al.l 
Il99d ). 

5.4 Obscuration 

A further step that we can take towards understanding the 
nature of the galaxies is to investigate the relative luminosi- 
ties at UV and submm wavelengths, which can give informa- 
tion about the fraction o f star-formation that is obscured by 
dust in the galaxies fe.g. lBuat et alj|2010l ; Iwiiesinghe et al.l 
l201lf ). The submm luminosity represents the energy ab- 
sorbed and re-radiated by dust. Naively, we might expect 
that this energy originated as UV radiation from young 
stars, hence the ratio of submm light to UV light detected 
is directly related to the fraction of UV light which is ob- 
scured by dust. If UV light is assumed to come primar- 
ily from star-forming regions, this is a measure of the ra- 
tio of obscured to unobscured star-formation. However both 
the UV and the submm radiation could also trace popula- 
tions unrelated to star-formation: there are open questions 
as to how much UV radiation can be produced by evolved 
stars (|Chavez fc Bertonel l201ll . and references therein) as 
well as how much of the dust probed at s ubmm wavelengths 
is heated by old stars in the galaxy dBendo et al.l l2010l : 
iBoselli et al]|2010al : lLaw. Gordon fc MisseltlbOllI). Detaile d 
radiative transfer calculations (e.g. Popescu et al.l |201lf ). 
which have been used to make detailed predictions for a 
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Figure 16. Stacked dust mass as a function of g — r colour, redshift and stellar mass. Dust mass is derived from equation J9]l using 
the fitted dust temperatures from Fig. 1141 Error bars include the statistical Icr errors in the bins as described in Section 13.31 with an 
additional contribution due to the error on the fitted temperature. The lensing contribution has been removed from the r ed bin s, and 
error bars include the associated uncertainty. The shaded region shows the range of characteristic dust masses measured bv lDlll . which 
evolve from 3.8 X 10 7 M at z < 0.1 to 2.2 X 10 s M atz~ 0.35. 
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Figure 17. Stacked dust mass per unit stellar mass as a function of g — r colour, redshift and stellar mass. Error bars include the 
statistical lcr errors in the bins as descri bed i n Section 13.31 with an additional contribution due to the error on the fitted temperature. 
The shaded region shows the range from ID 111 from 2 X 10 — 3 at z < 0.1 to 7 X 10~ 3 at z ~ 0.35. 



few well-studied spiral galaxies, can be used in the future to 
analyse statistical samples of galaxies to address this ques- 
tion in a quantitative way. For the moment, however, the 
generalisation to galaxy populations as a whole is uncertain. 

In Fig. [18] we stack L250/ Lnuv for the 7Vt/V-detected 
sample. Since we require NUV detections for this, we use 
the NUV — r colour which is likely to be a cleaner colour 
separation, and we use L250 instead of Ltir because the 
monochromatic luminosity is not model-dependent. Simu- 
lations showed that stacking this ratio is robust even for 
small 250 fim fluxes with low signal-to-noise, since the NUV 
fluxes all have reasonably high signal-to-noise (on the con- 



trary stacking Lnuv I L250 was found to be unreliable in 
simulations since this quantity diverges as the 250 ^m flux 
approaches zero). As before, we correct 250 /im fluxes for 
the expected contribution from lensing as described in Sec- 
tion 03] 

Some of the results implied by Fig. [18] are not trivial 
to explain, and should be treated with caution since there 
is a strong bias introduced by the UV selection. It appears 
that the obscuration increases with increasing stellar mass 
for blue galaxies, but there appears to be a decrease with 
redshift, at least for stellar masses < 5x 10 10 M . This con- 
trasts with the increase in 250 /jm luminosity with redshift, 
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which would imply that while the obscured SFR increases 
with redshift up to z — 0.3, the unobscured SFR (UV lumi- 
nosity) must increase faster for the relative obscuration to 
fall. However it is likely that these observations are affected 
by selection bias: we only detect the low mass galaxies in 
the UV if they are relatively unobscured, and as redshift 
increases we detect fewer and fewer of the obscured ones. 
This effect could cancel out any intrinsic increase in obscu- 
ration with redshift, and cause the observed L250/ 'Lnuv to 
decrease. 

In contrast we detect almost exactly the opposite trends 
in the red sample, and in the high mass end of the green 
sample, which suggests that the selection bias could hide 
similar trends in blue galaxies (which generally have much 
lower stellar mass for the same redshift). Number statistics 
are poor in the red bin because the selection is naturally 
biased against red galaxies, and errors are compounded by 
the uncertainty on the lensing contamination. Nevertheless 
the observed trends of increasing obscuration with increas- 
ing redshift and with decreasing stellar mass cannot be ex- 
plained by the selection bias. These trends are both per- 
fectly consistent with the trends in Md U8 t/M s tar in Fig. fTTl 
a higher amount of dust per stellar mass is almost certain to 
increase the obscuration of UV light. The red galaxy sample 
therefore appears to contain a larger fraction of obscured 
star-forming galaxies at higher reds hifts and lower m asses. 
This is consistent with the findings of lZhu et al.l (|201l|) in an 
analysis of the mid-IR colours of optically-selected galaxies 
at redshifts b e tween 0.1 and 0.5. It is also in agreement with 
iToieiro et al.l (|201ll 'l. who stacked SDSS spectra of luminous 
red galaxies (LRGs), and fitted stellar population models 
to obtain representative star-formation histories, metallicity 
and dust content as a function of colour, luminosity and 
redshift. Their results showed strong correlations of dust 
extinction with optical luminosity and redshift which are 
consistent with our own findings. Such agreement between 
independent measures of the dust extinction is encouraging. 

In any case we must be careful in the interpretation 
of L250/LNUV as a tracer of obscuration, in particular due 
to the potential for L250 to be uncorrelated with star for- 
mation. We showed in Section \5. II that the conversion from 
SPIRE luminosities to Ltir is extremely model-dependent, 
and to plot the ratio Ltir/Luv using our SED fits for 
Ltir would be misleading when Ltir is based only on 
the SPIRE photometry . In addition, it has been shown by 
IWiiesinghe et all (|201lh that the ratio of Ltir (from fit- 
ting SEDs to GAMA/H- ATLAS data including PACS and 
SPIRE) to Luv is poorly correlated with other measures of 
obscuration (the Balmer decrement and UV slope), proba- 
bly because they trace a different component of the dust in 
galaxies. We therefore hesitate to take this particular anal- 
ysis any further without the addition of shorter wavelength 
data to better constrain the IR SED. 



6 CONCLUSIONS 

We have conducted the first submm stacking analysis of a 
large sample of about 80,000 galaxies uniformly selected by 
optical (r-band) magnitude. We divided the sample by rest- 
frame colour, absolute magnitude/stellar mass and redshift 
(0.01 ^ 2 0.35) and stacked into SPIRE maps cover- 



ing about 126 square degrees at 250, 350 and 500 fim. We 
used a simple (but effective) deblending method to avoid the 
problem of over-estimating the flux of blended sources when 
stacking in confused images; this ensures that stacked flux 
ratios are not biased by the increasing level of confusion at 
longer wavelengths. Our main results are summarised below: 

(i) The submm fluxes of all but the most massive 
optically-selected galaxies are below the 5a limits of the 
H-ATLAS, yet with the large sample size made possible by 
the coverage of H-ATLAS and GAMA we are able to probe 
more than an order of magnitude below these limits using 
stacking. 

(ii) We estimate that the total emission from optically- 
selected galaxies at r < 19.8 and z < 0.35 accounts for only 
5.0 ± 0.4 % of the cosmic infrared background at 250 /im. At 

2 < 0.28, where the sample is complete to below M* , this 
fraction is 4.2L0.3 %. Of this, roughly 60 per cent originates 
from blue galaxies, and 20 per cent each from the red and 
green bins of our sample. 

(iii) We derive the total k- and e-corrected luminosity 
density of the Universe at z = to be 3.9 ± 0.3 x 10 33 , 
1.5 ± 0.1 x 10 33 and 4.3 ± 0.3 x 10 32 WMpc~ 3 h 70 at 250, 
350 and 500 respectively (/170 = -Ho/70 kms -1 Mpc -1 ). 

(iv) We show that stacked fluxes of red galaxies can be sig- 
nificantly contaminated by the lensing of background SMGs. 
Using models for the lensing amplification distribution and 
observed lens number counts, we estimate that around 10, 
20, and 30 per cent (at 250, 350 and 500 /im respectively) 
of the stacked fluxes is likely to result from lensing. We cor- 
rect our stacked results for this contamination to red galaxy 
stacks making reasonable assumptions for the redshift dis- 
tribution of lensed flux, and include the uncertainty from 
the lens number counts. 

(v) We observe a strong dependence of submm luminos- 
ity on optical colour (g — r) and stellar mass or M r , with red 
galaxies being up to an order of magnitude less luminous 
than blue galaxies of equal stellar mass. The luminosities of 
green galaxies are intermediate between the two. The ob- 
served trends of SPIRE luminosities are not strongly depen- 
dent on the SED model assumed, and cannot be explained 
by lensing, which implies a fundamental difference between 
the dust emission properties of red and blue galaxies. 

(vi) We measure cold dust temperatures that vary 
strongly as a function of stellar mass in blue galaxies, from 
~ 11 K at 3 x 10 8 M Q to ~ 28 K at 5 x 10 10 M Q . Correct- 
ing for the contamination from lensing, red galaxies have 
dust temperatures ~ 16 K at all stellar masses between 

3 x 10 9 - 8 x 10 10 M at 2 < 0.35. The dust temperatures of 
green galaxies appear to have a greater scatter (with mean 
T = 19.4 K) but are not correlated with stellar mass as with 
the blue; this is indicative of a mixed population. Temper- 
ature values depend on the assumption of a constant emis- 
sivity parameter ft = 2. 

(vii) The temperature variation can account for much of 
the difference in luminosities between red and blue galaxies; 
however it is not responsible for an increase in luminosity 
with redshift by a factor around 2 for blue galaxies at a 
given stellar mass. This appears to be due to an increase 
in the dust masses of galaxies of all stellar masses by a fac- 
tor of 3 — 4 between 2 ~ and 2 ~ 0.3. The red sample 
exhibit a stronger luminosity evolution, for which the likely 



© 0000 RAS, MNRAS 000, 000-000 



H- ATLAS/ GAM A: dust in optically selected galaxies 29 



BLUE: LOOS NUV-r<-0. 27-0.1 7M r GREEN: -0.27-0. 1 7M r 5 NUV-r< 1 .23-0. 1 7M r 



RED: 1 .23-0.1 7M r S NUV-r<7.00 




9 10 11 

l°9,o(M stor [MJ) 



9 10 1 

l°9,o(M stor [MJ) 



9 10 1 

l°g 10 (M slor [MJ) 



Figure 18. Stacked 250 fim/ NUV luminosity as a function of NUV — r colour, redshift and stellar mass. This fraction can be used 
as a proxy for the relative obscuration of star formation, subject to the limitations discussed in the text. Error bars are the statistical 
Ict errors in the bins as described in Section 13.31 Data and errors in the red bin incorporate the correction for lensing described in 
Section I3~2l 



explanation is also dust mass evolution. Due to the lensing 
uncertainty we cannot rule out evolution in the tempera- 
tures of the red sample, although there is no temperature 
evolution in the other colour bins. 

(viii) We fit the evolution in L250 (which is not dependent 
on the temperature) with the function L(z) oc (1 + z) a , and 
obtain indices for the three colour bins at Afg-tar — 2 — 7 x 
10 10 M . We find a = 4.0±0.2 for blue galaxies, a = 1.1±0.4 
for green and a =7.7±1.6 for red (the larger error on the 
evolution of red galaxies is due to the uncertainty on the 
lensing correction). Consistent rates of evolution are also 
derived for the dust masses. The evolution suggests a change 
in the dominant population of red galaxies, from passive 
systems at low redshift to obscured star-forming systems at 
higher redshift. 

(ix) The redshift evolution of galaxies classified as green 
seems to indicate a change in the nature of galaxies selected 
in this way, from a sample dominated by blue-cloud-like 
galaxies at low redshift and low M r (or M star ) to a sam- 
ple more similar to the red bin at the higher redshift and 
brighter M r . 

(x) Deriving TIR luminosities is problematic with only 
the SPIRE data, and we show that the results obtained de- 
pend sensitively on the SED model used (therefore the dust 
temperature) . The low temperatures implied by the SPIRE 
colours i ndicate that a cold model such as the H-ATLAS 
SED fits jSmith et alj|2011bf ) are more appropriate than ear- 
lier models based on IRAS and SCUBA data (CE01). 

(xi) The dust-to-stellar mass ratio is strongly anti- 
correlated with stellar mass, varying by more than an order 
of magnitude between M star ~ 10 s — 10 11 Mq. This relation- 
ship appears to vary little between different optical colours, 
although it evolves toward higher values with increasing red- 
shift. These results provide a challenge to dust formation 
models that rely on a purely stellar source of dust, implying 
a need for dust formation in supernovae and/or the ISM to 



reach the high dust masses in galaxies at the lower end of 
the stellar mass function. 

(xii) We attempt to explore the obscuration of galaxies 
in our sample using the L250 / 'Lmuv ratio, and see that red 
and green galaxies may become more obscured at increas- 
ing redshift and decreasing stellar mass (results for the blue 
galaxies are unclear due to selection bias). This conclusion is 
dependent on the assumption that this ratio is a good tracer 
of obscuration, but due to uncertainties in the heating mech- 
anism for cold dust this may not be valid. Nevertheless, such 
trends in obscuration are consistent with the trends of lu- 
minosity and dust /stellar mass. 

This study is the first of its kind and provides some 
tantalising glimpses of the characteristics of emission from 
dust in normal galaxies. Our understanding of the IR SED 
of optically-selected galaxies and of the obscuration of star 
formation would be greatly improved by the availability of 
data covering the peak of the SED. In a future study we 
hope to stack data from PACS and WISE in order to make 
a much more detailed analysis of the full SED. 
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APPENDIX A: ESTIMATING FLUXES OF 
BLENDED SOURCES BELOW THE NOISE 
LEVEL 

Al Deblending individual sources 

When stacking sources we must be careful not to over- 
count the flux in blended sources, which would lead to over- 
estimation of stacked fluxes especially in bins whose galax- 
ies are more clustered, and in the longer wavelength im- 
ages which have lower resolution. Since many sources are 
close to or below the noise level in the images it is impos- 
sible to model them from the images themselves, and since 
submm flux is poorly correlated with optical flux, we have no 
other prior information to base models on. We must there- 
fore make some simplifying assumptions in order to avoid 
over-counting. 

Consider that two or more sources may be blended, but 
we do not know the true flux of either or their brightness 
ratio. How can we decide how much of the blended flux 
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to attribute to each source? We first make the assumption 
that all sources are unresolved and therefore have a shape 
given by the PSF. We treat the image pixel-by-pixel and 
assume that the fractional contribution to a pixel from each 
of the nearby blended sources is dependent only on the dis- 
tance of that pixel from the source in question. The con- 
cept is visualised in one dimension in Fig. I All Panel (a) 
shows two sources A and B at positions xa and xb, whose 
fluxes are distributed across the image as }a(x) and /a (a;) 
(Jy pixel" 1 ). The sources are blended in the image and the 
measured data (solid line) is given by 

ftot(x) = f A (x) + f B (x) (Al) 

In Fig. lAlT b) the sources are modeled by PSFs of equal 
height (pa, Pb', thin black lines). We can convolve the image 
data with the PSF of A to give the function 

F p , A (x')=p A *ftot- (A2) 

To compute the convolution the PSF is shifted so that the 
peak is at the origin. The convolution is a function of the 
offset x' , and the source flux is strictly given by the value at 

X = XA, 

f P ,A = F p ,a(x a )/Xpa (A3) 

(and similarly for B). The division by the sum of the PSF 
squared normalises the flux. However, both f p ,A and f Pt B 
now contain too much flux because both include all of the 
flux that is blended. This blended flux would therefore be 
counted twice if the sources were stacked. Instead of doing 
this, the PSF pa can be weighted by the function 

g A {x) = pa(x)/[ PA (x) + p B (x)] (A4) 

which is simply the fractional contribution from the PSF pa 
at position x to the total pa + Pb- Thus we can replace pa 
and pb with 'deblended' PSFs 

q A (x) = g A (x)p A (x) (A5) 

(and similarly qs(x)) which are given by the thick grey lines 
in Fig. lAlf b). Convolving the image data with each of these 
deblended PSFs gives a more conservative estimate of the 
total flux: 

F qiA (x') = q A * /tot (A6) 
Again the flux is given by the value of the convolution at 

x' = X A - 

f q ,A=F gtA (XA)/X P 2 A (A7) 

The total (deblended) flux f q ,A+f q ,B is the same as the total 
input flux under the functions in Fig. lAlf a - ). whereas the 
total of fp.A + f P ,B is greater because blended flux has been 
double-counted. This deblending method always conserves 
total flux, whatever the ratio of the fluxes. 

On the other hand, the individual fluxes measured using 
equation (IA7[) are not exactly correct because blended flux 
is shared evenly between the two sources, whereas ideally 
it should be distributed according to the flux ratio of the 
sources. Hence in this example the recovered flux of A is too 
low and that of B too high (this effect is worsened by closer 
proximity of the sources). However with no prior information 
on the true flux ratio this is the best estimate that can be 
made. 




image position, x 

Figure Al. (a) Two simulated point sources A and B in a one- 
dimensional image, modeled with the same PSF but different nor- 
malisations, represented by the dashed and dash-dotted curves 
respectively. The total flux in the image as a function of position 
x is given by equation (I AH and is represented by the solid line, 
(b) The thin black lines are the PSFs, pa and pb, centred at 
x = 200 and x = 300 respectively. Both PSFs have width a = 50. 
The thick grey lines are the PSFs weighted for deblending, q A 
and qs, given by equation | |A5| I. (c) The reconstructed sources 
given by the image data (solid line in panel a) weighted by the 
PSFs in the middle panel. The thin black lines arc obtained using 
the unweighted PSFs p in equation | |A2| |. and the thick grey lines 
using the weighted PSFs q in equation (I A6 1 ■ 



To generalise this method to a two-dimensional image 
with multiply-blended sources we start with an image ar- 
ray of the same dimensions as the data image, filled with 
values of zero. To this we add a PRF for every source in 
the input catalogue, centred on the pixel where the source 
is located and interpolated from the PSF with a small off- 
set to correctly account for sub-pixel-scale positioning. In 
the region of an isolated source this image will be identical 
to the individual PRF, but where sources are blended it is 
equal to the sum of the PRFs (analogous to the sum of the 
thin black lines in Fig. lAlf b)). Thus all multiple blends are 
automatically accounted for, and the image we have con- 
structed is analogous to the denominator in equation (| A4f) 
(i.e. pa +Pb + ■•■)• For each source i we derive the weighting 
function gi(x,y) as the ratio of the PRF to a cutout region 
of our all-PRFs image, as in equation ()A4|) . In other words, 
the weight given to the flux in a pixel (x, y) is the value of 
the PRF of the target at (x,y) divided by the sum total of 
the contributions of all PRFs in that pixel. We measure the 
flux of each source by convolving the data image ( Jy pixel - ) 
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with the weighted PRF, as in equation (|A6|) . We tested the 
method in simulated maps with realistic source densities and 
using the PRFs of the three SPIRE bands. We found that 
the correct mean and median fluxes were always recovered 
when stacking, and that convolving with the PRF without 
any deblending always led to an overestimate of the median 
and mean fluxes. 

The deblending technique is carried out before any bin- 
ning, so all catalogue sources in the field are automatically 
deblended, not just those in the same bin as the target 
in question. We note that the method is essentially very 
similar to the 'global deblen ding' technique described by 
iKurczvnski fe Gawiser] |2010l ). which was demonstrated in 
that paper to minimise bias and variance in the stacks. 

A2 Comparison to a statistical approach 

The problem of stacking into confused maps is not a new 
one, and other methods for removing the excess flux due to 
blending have been used in the literature. The advantages of 
the method described above are that it automatically takes 
into account blending between objects in different bins, and 
also allows for the possibility of different bins having differ- 
ent amounts of blending (e.g. due to t he stronger cluster - 
ing of more massive and red galaxies - IZehavi et al.ll201ll ). 
To check that the results of our method are reasonable, we 
compare them to a simple statistical approach to calculate 
the average fraction of the stacked flu x which results from 
multiple counting o f blended sources (|Serieant et alJfeoOSt 
I Bourne et al .1120 111 ). Even a randomly distributed sample 
would lead to some excess signal from the superposition of 
multiple targets, but the average excess can be measured 
by simply stacking random positions in the map and sub- 
tracting this signal from the target stack. We do this already 
when we perform background subtraction (see Section f3. If) . 
If however there is any clustering in the sample, then the 
probability of a target being superimposed on another tar- 
get is greater than that of a random position being super- 
imposed on a target, so the excess signal from blending in 
the target stack is not f ully cancelled out by the background 
subtraction. Following ISerieant et al.l ||2008h . the fractional 
flux contribution due to clustering is given by 



F 



w(0)e 



-9 2 /2a 



2n9&9 



(A8) 



where n is the background source density, w(6) is the two- 
point angular correlation function of the sample and a is the 
width of the Gaussian beam profile in the map that we stack 
into. The function w(6) describes the excess probability of a 
background source appearing at an angular distance 9 from 
a target position, compar ed with a random distri bution. We 
measured this using the lLandv fc Szalavl (|l993l ) estimator 
which counts pairs between the data (D) and random (R) 
positions as a function of separation 9: 

, n . DD -2DR- RR . . . 

= ^ (A9) 

We counted pairs in 40 radial bins logarithmically spaced 
between 4 and 180 arcsec. The results, averaged over the 
three fields, are shown in Fig. IA2I together with a power- 
law fit by linear regression, described by w(9) = (0.012 ± 
0.001)6< (_o ' 76±0 ' O3) . This fit is in good agreement with prev i- 
ous results from SDSS r-limited data (|Connollv et al.|[2002f ). 



1.00 




Figure A2. The two-point angular correlation function of the 
GAMA catalogue used in this work, averaged over the G09, 
G12 and G15 fields. Error bars are the standard errors be- 
tween the values obtained in the three fields. A power-law fit 
by linear regression (with free slope and normalisation) gives 
uj(0) = (0.012 ± 0.001)6>(-°- 76±0 - 03 > . 



Table Al. Flux correction factors C = 1/(1 + F) based on equa- 
tion lA8t , compared with the typical / average ratios of deblended 
to non-deblended flux using equations l|A7| l and (jA3j . 





Statistical 


Mean (median) 




correction (C) 


deblending correction 


250 A*m 


0.9702 ± 0.0001 


0.937 (1.000) 


350 fj,m 


0.9550 ± 0.0002 


0.932 (0.987) 


500 fj,m 


0.9326 ± 0.0003 


0.989 (0.902) 



By substituting into equation (|A8[) the fit to w(9) and 
the beam sizes at 250, 350 and 500 /J,m (18, 25, and 35 arcsec 
FWHM), we obtained the fractional contribution to stacked 
fluxes due to blending. The corrected flux is then obtained 
by multiplying the stacked flux by C — 1/(1+F). Results are 
summarised in Table [ATI We see that this average statistical 
correction factor is broadly similar to the typical correction 
to individual fluxes using the deblending technique (equa- 
tion (IA7|) ). The deblending technique has the advantage of 
correcting the flux of each target individually, so however 
the sample is binned the appropriate average correction is 
made. This is not the case using the correlation function of 
the entire catalogue, since that only gives a single correction 
factor. It would be possible to calculate separate correction 
factors for each bin, using the cross-correlation be tween the 
bin and the full sample (as in lBourne et alj|201ll ). but the 
uncertainties would be significantly increased since each bin 
contains a relatively small number of objects. 



APPENDIX B: SIMULATIONS OF BIAS IN 
THE MEDIAN 

We simulated power-law distributions of fluxes with Gaus- 
sian noise added, and found that in certain cases the median 
measured flux (true flux plus noise) was biased high with 
respect to the median true flux, as a result of noise in the 
measured values. The amount of bias depends on (i) the flux 
limits of the distribution; (ii) the la noise level in relation to 
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the flux limits; and (iii) the slope of the power law describing 
the underlying flux distribution. The bias only becomes ap- 
parent when considering distributions with a median signal- 
to-noise less than 5a. 

In order to ascertain the level of bias that could be 
present in our stacks we must consider the shape of the 
underlying (true) flux distribution of sources in the stacks. 
In order to do so we must look at the distribution of fluxes 
in the much brighter H-ATLAS detected sample, which are 
not dominated by noise. The situation is helped considerably 
by the fact that we bin the sample according to M r and 
redshift, meaning that each bin is likely to have a limited 
range of fluxes with a distribution determined by the LF. 
To estimate the flux limits in a given bin we can look at the 
subm m fluxes of the galaxies with the highest optical fluxes. 
iDllI show the distribution of r magnitude and £250 in the 
SDP ID catalogue. At r < 16 the catalogue contains the 
full range of submm fluxes, which at a given r spans 1.3 dex 
in iS25o- We inspected the Phase 1 reliable IDs with good 
spectroscopic redshifts (250 (im sources matched to SDSS 
data using the same method as ISmith et all l2011al ; these 
will be described by Hoyos, C. et al. in preparation). We 
found the same range of 1.3 dex in r < 16 sources in Phase 
1 as in SDP. 

We therefore assume that any bin of M r and redshift 
will have fluxes within a limited range. The actual limits 
of this will depend on the range of 

M r and z sampled by the bin (although the majority of 
fluxes at all M r and z < 0.35 will lie within the range 
0.1-100 mjy). However, since the range in a bin is deter- 
mined by the LF, we can safely assume that the logarithmic 
range R = log 10 (S , ma x/Smin) = 1.3 will be consistent in all 
bins. Similar ranges of 1.3 dex are expected at all three 
SPIRE wavelengths (although of course fluxes at different 
wavelengths will be offset with respect to each other due to 
the shape of the SED). Within these limits the fluxes are 
assumed to follow a power-law distribution: in the Phase 1 
IDs this is approximately described by differential number 
counts dN/dS oc S -2 ' 5 , although this is unreliable due to the 
incompleteness of the ID catalogue. A similar slope is appar- 
ent in the 250 ji m number counts f rom P(D) analysis of the 
HerMES maps jGIenn et alj|20ich . although there is some 
evidence in that data for a shallower slope at S250 S; 10 mjy. 
We must however remember that the number counts in our 
bins will not follow the overall submm number counts, since 
our bins contain only a narrow distribution of M r and partic- 
ularly of redshift. In sufficiently narrow redshift bins the dis- 
tribution of the counts will approach the submm LF: this has 
a slope of -1.01 at t he faint end (L < V) of the H-ATLAS 
LFs derived bv lDlll . We can therefore be confident that in 
finite redshift bins the s l ope w ill be intermediate between 
— 1 and —2.5. lLapi et all (|201lh have modelled the number 
counts to fit the data from H-ATLAS, HerMES and BLAST. 
We split these into redshift bins (Sz = 0.05) and found that 
the faint-end differential counts at z < 0.35 follow a slope 
of approximately —2. Our redshift bins have similar widths 
(0.05 < Sz < 0.11) hence it is reasonable to assume that the 
flux distributions in the bins will have a similar slope. 

For these reasons we conclude that a simulated flux 
distribution that spans a range R — 1.3 dex with a power 
law dN/dS oc S~ 2 is representative of our bins. We there- 
fore created 16 simulated distributions with these param- 



eters, with a range of lower limits between 0.1 — 30mJy: 
the corresponding upper limits are 20 times larger, given 
by R = 1.3. The resulting distributions have medians in the 
range 0.3 — 70 mjy, which is sufficient to cover all the bins in 
our real data. We added random noise to these true flux dis- 
tributions, where the noise values were drawn from a Gaus- 
sian distribution with zero mean, and a given by the average 
total noise level (instrumental plus confusion) in the Phase 
1 maps: a = 6.7, 7.9, 8.8 mjy beam" 1 at 250, 350, 500 
respectively. We then compared median measured flux (true 
flux plus noise) to the median true flux, and quantified the 
bias factor as the ratio of measured to true median flux. 

We corrected the measured median fluxes in our stacked 
data by interpolating the relationship between true median 
flux and measured median flux from the simulations. This 
was done separately for each of the three bands (since the 
bias behaves differently as a result of the different noise lev- 
els). All correction factors are in the range 0.6 — 1.0, which is 
generally small in comparison to the range of stacked fluxes 
resulting from true differences between the bins. 

For these corrections we have assumed that the slope of 
differential number counts is —2 and that fluxes in each bin 
lie in a range given by R = 1.3 dex. However, the bias factor 
depends on both of these variables, as we show in Fig. HU 
To test the sensitivity of our results to these parameters, we 
tried correcting our results by the bias factors obtained using 
different values. We tried slopes of dN/dS ranging from —0.5 
to —4 and R- values between 1 and 2. The level of bias varies, 
indicating that there is some uncertainty in the correction, 
but we note that the corrections for slopes between —1.5 
and —2 give almost identical corrections (for R = 1.3). We 
are sure that the slope is between —1 and —2.5, and over 
this range the bias corrections do not vary by more than 
15%. Increasing the range R has a greater impact on the 
correction; however the value of 1.3 is well motivated and 
in narrow bins of M r and redshift there is good reason to 
expect the flux range to be limited. 

Crucially, all of our conclusions remain valid, and all 
trends remain significant, for any combination of slope (—0.5 
to —4) and range (1 to 2). This is equally true if we make 
no corrections. 

An alternative to these corrections would be to use the 
mean instead of the median when stacking. The mean is 
not altered by the effects of noise as we found the median 
to be; however the mean will be highly biased simply by 
the skewed shape of the distribution. In fact the bias in the 
mean at all flux levels is equal to the maximum bias seen 
in the median at the lowest flux levels. The reason for this 
is simple: at the lowest flux levels, where noise dominates 
over the true flux, the distribution of measured fluxes closely 
resembles the noise distribution - a symmetrical Gaussian - 
but instead of being centred on zero as the noise is, it has the 
same mean as the true flux distribution. The mean is not 
altered by the addition of noise if the mean noise is zero. 
The measured flux distribution is therefore symmetrical in 
this case, hence the median is equal to the mean. Thus the 
maximum amount of bias in the median occurs when the 
shape of the distribution becomes dominated by the noise 
rather than by the true fluxes. 

In conclusion, we choose to use the corrected median 
estimator rather than the mean, because at all but the low- 
est fluxes the median is a better descriptor of the underlying 
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(true) flux distribution. At the lowest fluxes the median be- 
comes biased and approaches the mean of the distribution. 
However as a result of our binning scheme we can make a 
good estimate of the shape of the underlying flux distribu- 
tion and can be reasonably certain of the correction factors 
for the bias. The fact that all of our ultimate conclusions 
remain valid for any reasonable choice of the distribution 
indicates the robustness of our results to these corrections. 



APPENDIX C: POSTAGE STAMPS, 
HISTOGRAMS AND SPECTRAL ENERGY 
DISTRIBUTIONS OF STACKS 

In Figures [C1I - IC4I we choose some example bins to show the 
stacked postage-stamp images, the distribution of measured 
fluxes, and the SED data and model. The Figures below 
contain the following information: 

Postage stamps of the stack in the three SPIRE bands 
are shown with contours at signal-to-noise levels of 5, 10, 
15, 20, 25, 50, 100, 150, 200 & 250. The images are each 41 
pixels square, corresponding to 3'25" at 250 ^tm and 6'50" 
in the other two bands. These images are illustrative only, 
and were made by stacking in the PSF-filtered, background- 
subtracted SPIRE maps (we do not measure fluxes in these 
maps but using the deblend filter method described in Ap- 
pendix The flux and signal-to-noise reached in the cen- 
tral pixel of each postage stamp agree with the stacked val- 
ues in Table [TJ although the postage stamps show slightly 
boosted fluxes due to blending not being accounted for. Sim- 
ilar agreement was found in all stacks, although only a se- 
lection are shown here for brevity. 

We also show histograms of the measured fluxes in the 
stack (red) and of a set of fluxes measured at random posi- 
tions in the background (blue), as described in Section [4] • 
The number shown is the KS probability that these two sam- 
ples were drawn from the same distribution. Beneath these is 
plotted the single-component SED fitted to the three stacked 
fluxes (plotted in the rest frame), assuming (8 = 2. The SED 
is fitted to obtain the temperature which is printed over the 
SED, as described in Section [5] 



© 0000 RAS, MNRAS 000, 000-000 




Figure Bl. Results of a set of stacking simulations: (a) Comparing median measured flux and median true flux, for distributions in 
various flux ranges. The range in each simulation is given by R = log 1 Q(S' max /S m i n ) = 1.3 (the length of the grey horizontal bar), and the 
simulations were run 16 times with different S m i n , S maJt values to produce distributions with a range of median (true) fluxes. Each coloured 
line in the Figure connects 16 data points, one point from each simulation, showing how the median measured flux (after adding noise) 
varies as the median true flux is varied. Noise is drawn from a Gaussian distribution centred on zero, with a = 6.7, 7.9, 8.8 mjy beam -1 
at 250, 350, 500 ^m respectively. Different coloured lines correspond to sets of simulations with different slopes of dN/dS between — 1 
and —4. (b) Bias factor (ratio of median measured to median true flux) as a function of median true flux, for the same set of simulations, 
showing how the bias depends on the slope of dN/dS for R = 1.3. Note that bias does not vary monotonically with slope for fixed R, but 
is greatest for a slope of —1.5 in this case, (c) Bias factor as a function of median true flux, showing the effect of increasing the range 
of true fluxes. Different lines correspond to sets of simulations with different values of R between 1 and 2; in each of these the slope of 
dN/dS is -2. 
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Figure CI. A stack of galaxies with blue g — r colours, median 
M r = —21.1, median z = 0.11. This stack is one of the brightest 
in submm flux, and contains 1567 objects. Stacked, PSF-filtered 
postage stamps are shown with a logarithmic greyscale between 
0.0001 — 0.05 Jy beam -1 and signal-to-noise contours at 5, 10, 15, 
20, 25, 50, 100, 150, 200, 250. Histograms of the measured fluxes 
in the stack (red) and in a random background stack (blue) are 
shown with the KS probability that they were drawn from the 
same distribution. Also plotted is the rest-frame SED fit with the 
temperature indicated, assuming = 2. More details are given in 
the text of Appendix [C] 
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Figure C2. A stack of galaxies with green g — r colours, median 
M r = —21.1, median z = 0.15. This stack has moderate submm 
flux, and contains 570 objects. Plots arc as in Fig. IC1I except 
that the postage stamps are plotted on a logarithmic greyscale 
between 0.0001 - 0.02 Jy beam -1 . 
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Figure C3. A stack of galaxies with red g — r colours, me- 
dian M r = —21.7, median z = 0.21. This stack has among the 
faintest submm fluxes, and contains 1111 objects. Plots are as in 
Fig. ICll cxccpt that the postage stamps arc plotted on a logarith- 
mic greyscale between 0.0001 — 0.005 Jy beam -1 . 
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Figure C4. A stack of galaxies with red g — r colours, me- 
dian M r = —22.5, median z = 0.27. This stack has among the 
faintest submm fluxes, and contains 1002 objects. Plots are as in 
Fig. ICll cxccpt that the postage stamps arc plotted on a logarith- 
mic greyscale between 0.0001 — 0.004 Jy beam -1 . 
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